#ifndef __MATH_PPAL_AXIS__
#define __MATH_PPAL_AXIS__

#include "../geometry/geometry_types.h"
#include "math_common.h"

#include <cassert>

//_________________________________________________________

    /**
     * PrincipalAxes3d enables the center and inertia axes of
     * a cloud of 3d points to be computed.
     */

	/* usage: 

	PrincipalAxes3d axes ;
	axes.begin_points() ;
	FOR_EACH_VERTEX(MapComponent, comp, it) {
		const Point3f& p = it->point() ;
		axes.point(p) ;
	}
	axes.end_points() ;

	Point3f center = axes.center() ;
	Vector3f Z0 = axes.axis(0) ;
	Vector3f Z1 = axes.axis(1) ;
	Vector3f Z2 = axes.axis(2) ;

	*/


    class MATH_API PrincipalAxes3d 
	{
    public:
        PrincipalAxes3d() ;
        void begin_points() ;
        void end_points() ;
        void add_point(const Point3f& p, double weight = 1.0) ;
        
        Point3f center() const ;
        const Vector3f& axis(int i) const ;
        double eigen_value(int i) const ; 
        
        
    private:
        double		center_[3] ;
        Vector3f	axis_[3] ;
        double		eigen_value_[3] ;
        
        double		M_[6] ;
        int			nb_points_ ;
        double		sum_weights_ ;
    } ;
    
//_________________________________________________________

    /**
     * PrincipalAxes2d enables the center and inertia axes of
     * a cloud of 2d points to be computed.
     */
    class MATH_API PrincipalAxes2d {
    public:
        PrincipalAxes2d() ;
        void begin_points() ;
        void end_points() ;
        void add_point(const Point2f& p, double weight = 1.0) ;
        
        Point2f center() const ;
        const Vector2f& axis(int i) const ;
        double eigen_value(int i) const ; 
        
        
    private:
        double		center_[2] ;
        Vector2f	axis_[2] ;
        double		eigen_value_[2] ;
        
        double		M_[3] ;
        int			nb_points_ ;
        double		sum_weights_ ;
    } ;

//_________________________________________________________

    inline Point3f PrincipalAxes3d::center() const {
        return Point3f(center_[0], center_[1], center_[2]) ;
    }

    inline const Vector3f& PrincipalAxes3d::axis(int i) const {
        assert(i >= 0 && i < 3) ;
        return axis_[i] ;
    }

    inline double PrincipalAxes3d::eigen_value(int i) const {
        assert(i >= 0 && i < 3) ;
        return eigen_value_[i] ;
    }

    //______________________________________________________

    inline Point2f PrincipalAxes2d::center() const {
        return Point2f(center_[0], center_[1]) ;
    }

    inline const Vector2f& PrincipalAxes2d::axis(int i) const {
        assert(i >= 0 && i < 2) ;
        return axis_[i] ;
    }

    inline double PrincipalAxes2d::eigen_value(int i) const {
        assert(i >= 0 && i < 2) ;
        return eigen_value_[i] ;
    }

#endif

