/****************************************************************************
* VCGLib                                                            o o     *
* Visual and Computer Graphics Library                            o     o   *
*                                                                _   O  _   *
* Copyright(C) 2004                                                \/)\/    *
* Visual Computing Lab                                            /\/|      *
* ISTI - Italian National Research Council                           |      *
*                                                                    \      *
* All rights reserved.                                                      *
*                                                                           *
* This program is free software; you can redistribute it and/or modify      *   
* it under the terms of the GNU General Public License as published by      *
* the Free Software Foundation; either version 2 of the License, or         *
* (at your option) any later version.                                       *
*                                                                           *
* This program is distributed in the hope that it will be useful,           *
* but WITHOUT ANY WARRANTY; without even the implied warranty of            *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
* for more details.                                                         *
*                                                                           *
****************************************************************************/

#ifndef VCG_USE_EIGEN
#include "deprecated_point4.h"
#else

#ifndef __VCGLIB_POINT4
#define __VCGLIB_POINT4

#include "../math/eigen.h"

namespace vcg{
template<class Scalar> class Point4;
}

namespace Eigen{
template<typename Scalar>
struct ei_traits<vcg::Point4<Scalar> > : ei_traits<Eigen::Matrix<Scalar,4,1> > {};
}

namespace vcg {
/** \addtogroup space */
/*@{*/
    /**
        The templated class for representing a point in 4D space.
        The class is templated over the ScalarType class that is used to represent coordinates. 
				All the usual operator (* + - ...) are defined. 
     */

template <class T> class Point4 : public Eigen::Matrix<T,4,1>
{
	typedef Eigen::Matrix<T,4,1> _Base;
	using _Base::coeff;
	using _Base::coeffRef;
	using _Base::setZero;
	using _Base::data;
	using _Base::V;

public:

	_EIGEN_GENERIC_PUBLIC_INTERFACE(Point4,_Base);
	typedef Scalar ScalarType;

	VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point4)
	
	enum {Dimension = 4};

	inline Point4() : Base() {}
	inline Point4( const T nx, const T ny, const T nz , const T nw ) : Base(nx,ny,nz,nw) {}
	inline Point4(const T p[4]) : Base(p) {}
	inline Point4(const Point4& p) : Base(p) {}
	template<typename OtherDerived>
	inline Point4(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
	

	/// importer from different Point4 types
	template <class Q> inline void Import( const Point4<Q> & b ) { *this = b.template cast<T>(); }
	
	/// constuctor that imports from different Point4 types
  template <class Q>
  static inline Point4 Construct( const Point4<Q> & b ) { return b.template cast<T>(); }


//@{

	inline T &X() {return Base::x();}
	inline T &Y() {return Base::y();}
	inline T &Z() {return Base::z();}
	inline T &W() {return Base::w();}

	inline const T & V ( const int i ) const
	{
		assert(i>=0 && i<4);
		return data()[i];
	}
	inline T & V ( const int i )
	{
		assert(i>=0 && i<4);
		return data()[i];
	}
	
	/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
	/// Useful for managing in a consistent way object that could have point2 / point3 / point4
	inline T Ext( const int i ) const
	{
		if(i>=0 && i<=3) return data()[i];
		else             return 0;
	}
//@}
	
	inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const
	{	
		Point4 res;
		const Point4 &y = *this;

		res[0] =  y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+
		          x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2];
		res[1] =  y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]*
		          z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2];
		res[2] = -y[0]*z[1]*x[3]+x[0]*z[1]*y[3]+y[0]*x[1]*
  		        z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3];
		res[3] = -z[0]*y[1]*x[2]-y[0]*x[1]*z[2]+x[0]*y[1]*
		          z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2];
		return res;
	}
	
	/// Homogeneous normalization (division by W)
	inline Point4 & HomoNormalize() {
		if (data()[3]!=0.0) {	Base::template start<3>() /= coeff(3); coeffRef(3) = 1.0; }
		return *this;
	}
	
//@{
  /** @name Comparison operators (lexicographical order)
  **/
	inline bool operator <  ( Point4 const & p ) const
	{
		return	(data()[3]!=p.data()[3])?(data()[3]<p.data()[3]):
				(data()[2]!=p.data()[2])?(data()[2]<p.data()[2]):
				(data()[1]!=p.data()[1])?(data()[1]<p.data()[1]):
				(data()[0]<p.data()[0]);
	}
	inline bool operator >  ( const Point4 & p ) const
	{
		return	(data()[3]!=p.data()[3])?(data()[3]>p.data()[3]):
				(data()[2]!=p.data()[2])?(data()[2]>p.data()[2]):
				(data()[1]!=p.data()[1])?(data()[1]>p.data()[1]):
				(data()[0]>p.data()[0]);
	}
	inline bool operator <= ( const Point4 & p ) const
	{
		return	(data()[3]!=p.data()[3])?(data()[3]< p.data()[3]):
				(data()[2]!=p.data()[2])?(data()[2]< p.data()[2]):
				(data()[1]!=p.data()[1])?(data()[1]< p.data()[1]):
				(data()[0]<=p.data()[0]);
	}
	inline bool operator >= ( const Point4 & p ) const
	{
		return	(data()[3]!=p.data()[3])?(data()[3]> p.data()[3]):
				(data()[2]!=p.data()[2])?(data()[2]> p.data()[2]):
				(data()[1]!=p.data()[1])?(data()[1]> p.data()[1]):
				(data()[0]>=p.data()[0]);
	}
//@}
	
//@{
  /** @name Dot products
  **/

	inline Point4 operator ^ (  const Point4& p ) const
	{
		assert(0 && "not defined by two vectors (only put for metaprogramming)");
		return Point4();
	}

	/// slower version, more stable (double precision only)
	T StableDot ( const Point4<T> & p ) const
	{
		T k0=data()[0]*p.data()[0],	k1=data()[1]*p.data()[1], k2=data()[2]*p.data()[2], k3=data()[3]*p.data()[3];
		int exp0,exp1,exp2,exp3;

		frexp( double(k0), &exp0 );frexp( double(k1), &exp1 );
		frexp( double(k2), &exp2 );frexp( double(k3), &exp3 );

		if (exp0>exp1) { math::Swap(k0,k1); math::Swap(exp0,exp1); }
		if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); }
		if (exp0>exp2) { math::Swap(k0,k2); math::Swap(exp0,exp2); }
		if (exp1>exp3) { math::Swap(k1,k3); math::Swap(exp1,exp3); }
		if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); }

		return ( (k0 + k1) + k2 ) +k3;
	}  
//@}


}; // end class definition


	/// slower version of dot product, more stable (double precision only)
template<class T>
double StableDot ( Point4<T> const & p0, Point4<T> const & p1 )
{
	return p0.StableDot(p1);
}  

typedef Point4<short>  Point4s;
typedef Point4<int>	   Point4i;
typedef Point4<float>  Point4f;
typedef Point4<double> Point4d;

/*@}*/
} // end namespace
#endif

#endif