/*
								+----------------------------------+
								|                                  |
								|    ***  Polygon template  ***    |
								|                                  |
								|   Copyright © -tHE SWINe- 2007   |
								|                                  |
								|            Polygon.h             |
								|                                  |
								+----------------------------------+
*/

#pragma once
#ifndef __POLYGON_INCLUDED
#define __POLYGON_INCLUDED

/**
 *	@file lml/Polygon.h
 *	@brief a simple polygon template
 *	@date 2007
 *	@author -tHE SWINe-
 *
 *	@date 2007-02-28
 *
 *	removed error in CPolygon2<TVertStruct>::Cut and CPolygon2<TVertStruct>::Split
 *	in cases edge was longer than 1 treshold for adding new edge vertex was bigger
 *	than epsilon so potencial maximal vertex - plane distance deviation was greater
 *	than epsilon. old code also always generated new vertex even in case it was not
 *	going to be used.
 *	fixed unnecessary / double vertices issue in CPolygon2<TVertStruct>::Split
 *
 *
 *	@date 2013-02-18
 *
 *	Fixed issues with overlap test and split polygon by plane.
 *	Added documentation comments and converted existing ones to doxygen format.
 *
 */

#include "../NewFix.h"
#include "../CallStack.h"
#include "../Vector.h"
#include <vector>
#include <algorithm>
#include "../Unused.h"
#include "../StlUtils.h"
#include "../MinMax.h"

/**
 *	@brief template for creating reference vertices for any arbitrary vertices
 *	@tparam TBaseVertex is a base vertex type
 */
template <class TBaseVertex>
struct TRefVertex {
	TBaseVertex *m_p_ref; /**< @brief referenced vertex */
	std::vector<TBaseVertex*> &m_r_vertex_pool; /**< @brief vertex pool (here we going to add new vertices) */

	typedef TBaseVertex TIntermediate; /**< @brief intermediate type for calculations will be TBaseVertex */

	/**
	 *	@brief default constructor; creates a copy of the original vertex and adds it to the pool
	 *
	 *	@param[in] r_t_vertex is position of the vertex
	 *	@param[in] r_vertex_pool is vertex pool to store the vertices in
	 *
	 *	@note This function throws std::bad_alloc.
	 */
	inline TRefVertex(const TBaseVertex &r_t_vertex, std::vector<TBaseVertex*> &r_vertex_pool) // throw(std::bad_alloc)
		:m_p_ref(0), m_r_vertex_pool(r_vertex_pool)
	{
		m_p_ref = new TBaseVertex(r_t_vertex); // alloc a new TBaseVertex
		m_r_vertex_pool.push_back(m_p_ref); // insert it to the list
	}

	/**
	 *	@brief constructor (reference vertex that is already in the list)
	 *
	 *	@param[in] p_vertex_it is an iterator pointing to the vertex
	 *	@param[in] r_vertex_pool is vertex pool to store the vertices in
	 */
	inline TRefVertex(typename std::vector<TBaseVertex*>::iterator p_vertex_it,
		std::vector<TBaseVertex*> &r_vertex_pool)
		:m_p_ref(*p_vertex_it), m_r_vertex_pool(r_vertex_pool)
	{}

	/**
	 *	@brief copy constructor
	 *	@param[in] r_t_vertex is a reference vertex to be copied
	 */
	inline TRefVertex(const TRefVertex &r_t_vertex)
		:m_p_ref(r_t_vertex.m_p_ref), m_r_vertex_pool(r_t_vertex.m_r_vertex_pool)
	{}

	/**
	 *	@brief copy operator (need it to avoid unwanted conversions)
	 *	@param[in] r_t_vertex is a reference vertex to be copied
	 *	@return Returns a reference to this.
	 */
	inline TRefVertex &operator =(const TRefVertex &r_t_vertex)
	{
		m_p_ref = r_t_vertex.m_p_ref;
		m_r_vertex_pool = r_t_vertex.m_r_vertex_pool;

		return *this;
	}

	/**
	 *	@brief conversion from TBaseVertex to TRefVertex (solved by adding vertex to the list)
	 *	@param[in] r_t_vertex is a vertex coordinate to be copied
	 *	@return Returns a reference to this.
	 *	@note This function throws std::bad_alloc.
	 */
	inline TRefVertex &operator =(const TBaseVertex &r_t_vertex) // throw(std::bad_alloc)
	{
		m_p_ref = new TBaseVertex(r_t_vertex); // alloc a new TBaseVertex
		m_r_vertex_pool.push_back(m_p_ref); // insert it to the list

		return *this;
	}

	/**
	 *	@brief means of conversion to base vertex type
	 *	@return Returns the coordinates of this vertex.
	 */
	inline operator TBaseVertex() const
	{
		return *m_p_ref;
	}

	/**
	 *	@brief means of conversion to raw vertex position
	 *	@return Returns the position of this vertex.
	 */
	inline operator Vector3f() const
	{
		return m_p_ref->v_position;
	}

	/**
	 *	@brief assigns vertex position
	 *	@return Returns a reference to this.
	 */
	inline TRefVertex &operator =(const Vector3f &r_v_position)
	{
		m_p_ref->v_position = r_v_position;
		return *this;
	}

	/**
	 *	@brief compares vertex coordinates
	 *	@param[in] r_t_vertex is a vertex coordinate to compare to
	 *	@return Returns true if the vertices are close enough, otherwise returns false.
	 *	@note The f_epsilon constant is used as maximal coordinate distance tolerance.
	 */
	inline bool operator ==(const TBaseVertex &r_t_vertex) const
	{
		return *m_p_ref == r_t_vertex;
	}

	/**
	 *	@brief vertex scaling operator
	 *	@param[in] f_scalar is a scaling factor
	 *	@return Returns a TBaseVertex because new vertex originated
	 *		(in case it will need to be converted to TRefVertex, it will be added to the list).
	 */
	TBaseVertex operator *(float f_scalar) const
	{
		return *m_p_ref * f_scalar;
	}

	/**
	 *	@brief vertex addition operator
	 *	@param[in] r_vertex is vertex coordinate to add
	 *	@return Returns a TBaseVertex because new vertex originated
	 *		(in case it will need to be converted to TRefVertex, it will be added to the list).
	 */
	TBaseVertex operator +(const TBaseVertex &r_vertex) const
	{
		return *m_p_ref + r_vertex;
	}

	/**
	 *	@brief vertex scaling operator
	 *	@param[in] f_scalar is a scaling factor
	 *	@return Returns a reference to this (no new vertices added).
	 */
	TRefVertex &operator *=(float f_scalar)
	{
		*m_p_ref *= f_scalar;
		return *this;
	}

	/**
	 *	@brief vertex addition operator
	 *	@param[in] r_vertex is vertex coordinate to add
	 *	@return Returns a reference to this (no new vertices added).
	 */
	TRefVertex operator +=(const TBaseVertex &r_vertex)
	{
		*m_p_ref += r_vertex;
		return *this;
	}
};
// todo don't really like this, put this away, at least proper pool

/**
 *	@brief namespace for utility classes for the CPolygon2 template
 */
namespace poly2 {

/**
 *	@brief small function object, finds a vertex forming a concave edge
 *	@tparam _TyVector is a vector type (specialization of Vector2 or Vector3)
 */
template <class _TyVector>
class CFindConcave;

/**
 *	@brief small function object, finds a vertex forming a concave edge
 *		(specialization for 3D vertices)
 *	@tparam _TyScalar is a scalar type
 */
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER >= 1400
template <class _TyScalar>
class CFindConcave<Vector3<_TyScalar> > {
#else // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
template <>
class CFindConcave<Vector3f> { // MSVC6 doesn't support advanced template specialization
	typedef float _TyScalar;
#endif // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
protected:
	Vector3<_TyScalar> m_v_prev2; /**< @brief the second previous vertex */
	Vector3<_TyScalar> m_v_normal; /**< @brief polygon normal vector */
	Vector3<_TyScalar> m_v_u; /**< @brief temporary edge vector, reused between the iterations */
	_TyScalar m_f_epsilon_ex; /**< @brief maximal allowed concavity */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] v_prev is the previous vertex
	 *	@param[in] v_prev2 is the second previous vertex
	 *	@param[in] v_normal is polygon normal vector
	 *	@param[in] f_epsilon_ex is maximal allowed concavity
	 */
	inline CFindConcave(Vector3<_TyScalar> v_prev, Vector3<_TyScalar> v_prev2,
		Vector3<_TyScalar> v_normal, _TyScalar f_epsilon_ex)
		:/*m_v_prev(v_prev),*/ m_v_prev2(v_prev2),
		m_v_normal(v_normal), m_f_epsilon_ex(f_epsilon_ex)
	{
		m_v_u = v_prev - v_prev2;
		m_v_u.Normalize();
	}

	/**
	 *	@brief tests a next vertex
	 *	@param[in] r_v_vertex_pos is vertex, coming after the previous ones specified
	 *		in the constructor, or after the one at the last call to this function
	 *	@return Returns true if the vertex forms a concave edge, otherwise returns false.
	 */
	inline bool operator ()(const Vector3f &r_v_vertex_pos)
	{
		Vector3<_TyScalar> v_v = m_v_prev2 - r_v_vertex_pos;
		v_v.Normalize();

		if(m_v_u.v_Cross(m_v_normal).f_Dot(v_v - m_v_u) < -m_f_epsilon_ex)
			return true;
		// norm of cross product is 1 since normal is perp, and we are projecting the distance
		// of the next edge: it is distance then (the edge height), not the angle cosine

		m_v_prev2 = r_v_vertex_pos;
		m_v_u = v_v;
		// shift vectors

		return false;
	}
};

/**
 *	@brief small function object, finds a vertex forming a concave edge
 *		(specialization for 2D vertices)
 *	@tparam _TyScalar is a scalar type
 *	@note Not sure which direction it marks as concave, need to test it.
 */
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER >= 1400
template <class _TyScalar>
class CFindConcave<Vector2<_TyScalar> > {
#else // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
template <>
class CFindConcave<Vector2f> { // MSVC6 doesn't support advanced template specialization
	typedef float _TyScalar;
#endif // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
protected:
	Vector2<_TyScalar> m_v_prev2; /**< @brief the second previous vertex */
	Vector2<_TyScalar> m_v_u; /**< @brief temporary edge vector, reused between the iterations */
	_TyScalar m_f_epsilon_ex; /**< @brief maximal allowed concavity */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] v_prev is the previous vertex
	 *	@param[in] v_prev2 is the second previous vertex
	 *	@param[in] v_normal is polygon normal vector
	 *	@param[in] f_epsilon_ex is maximal allowed concavity
	 */
	inline CFindConcave(Vector2<_TyScalar> v_prev,
		Vector2<_TyScalar> v_prev2, _TyScalar f_epsilon_ex)
		:m_v_prev2(v_prev2), m_f_epsilon_ex(f_epsilon_ex)
	{
		m_v_u = v_prev - v_prev2;
		m_v_u.Normalize();
	}

	/**
	 *	@brief tests a next vertex
	 *	@param[in] r_v_vertex_pos is vertex, coming after the previous ones specified
	 *		in the constructor, or after the one at the last call to this function
	 *	@return Returns true if the vertex forms a concave edge, otherwise returns false.
	 */
	inline bool operator ()(const Vector2f &r_v_vertex_pos)
	{
		Vector2<_TyScalar> v_v = m_v_prev2 - r_v_vertex_pos;
		v_v.Normalize();

		if(m_v_u.v_Orthogonal().f_Dot(v_v - m_v_u) < -m_f_epsilon_ex)
			return true;
		// norm of orthogonal vector in 2D is 1, and we are projecting the distance
		// of the next edge: it is distance then (the edge height), not the angle cosine

		m_v_prev2 = r_v_vertex_pos;
		m_v_u = v_v;
		// shift vectors

		return false;
	}
};

/**
 *	@brief simple function object for finding maximal
 *		absolute vertex deviation from normal plane
 *	@tparam _TyVector is a vector type (specialization of Vector2 or Vector3)
 */
template <class _TyVector>
class CFindMaxPlaneDeviation;

/**
 *	@brief simple function object for finding maximal absolute vertex
 *		deviation from normal plane (specialization for 3D vertices)
 *	@tparam _TyScalar is a scalar type
 */
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER >= 1400
template <class _TyScalar>
class CFindMaxPlaneDeviation<Vector3<_TyScalar> > {
#else // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
class CFindMaxPlaneDeviation<Vector3f> { // MSVC6 doesn't support advanced template specialization
	typedef float _TyScalar;
#endif // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
protected:
	float m_f_max_deviation; /**< @brief maximal deviation recorded */
	const Plane3<_TyScalar> &m_r_plane; /**< @brief polygon normal plane */

public:
	/**
	 *	@brief default constructor
	 *	@param[in] r_plane is polygon normal plane
	 */
	inline CFindMaxPlaneDeviation(const Plane3<_TyScalar> &r_plane)
		:m_r_plane(r_plane)
	{}

	/**
	 *	@brief updates maximal absolute vertex deviation from normal plane
	 *	@param[in] r_v_vertex_pos is a polygon vertex
	 */
	inline void operator ()(const Vector3<_TyScalar> &r_v_vertex_pos)
	{
		_TyScalar f_deviation = _TyScalar(fabs(m_r_plane.f_Vector_Dist(r_v_vertex_pos)));
		if(m_f_max_deviation < f_deviation)
			m_f_max_deviation = f_deviation;
	}

	/**
	 *	@brief gets the maximal absolute vertex deviation
	 *	@return Returns the maximal absolute vertex deviation.
	 */
	inline operator _TyScalar() const
	{
		return m_f_max_deviation;
	}
};

/**
 *	@brief small function object used to find third edge longer than epsilon
 *	@tparam _TyVector is a vector type (specialization of Vector2 or Vector3)
 */
template <class _TyVector>
class CFindLongEdges {
protected:
	typedef typename _TyVector::_TyScalar _TyScalar; /**< @brief a scalar type */
	_TyVector m_v_prev; /**< @brief the previous vertex */
	_TyScalar m_f_epsilon; /**< @brief long edge threshold */
	int m_n_edge_num; /**< @brief long edge counter @note The use of int is sufficient, only counts to two. */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] v_prev is the previous vertex
	 *	@param[in] f_epsilon_ex is long edge threshold
	 */
	inline CFindLongEdges(_TyVector v_prev, _TyScalar f_epsilon_ex)
		:m_v_prev(v_prev), m_f_epsilon(f_epsilon_ex), m_n_edge_num(0)
	{}

	/**
	 *	@brief finds long edges
	 *	@param[in] v_vertex is vertex, coming after the previous one specified
	 *		in the constructor, or after the one at the last call to this function
	 *	@return Returns true if it encountered at least three long edges, otherwise returns false.
	 */
	inline bool operator ()(_TyVector v_vertex)
	{
		if((m_v_prev - v_vertex).f_Length() >= m_f_epsilon) {
			if(m_n_edge_num == 2)
				return true;
			// in case there are three edges, triangle is not tiny

			++ m_n_edge_num;
		}
		// is the edge longer than epsilon?

		m_v_prev = v_vertex;

		return false;
	}
};

/**
 *	@brief small function object for evaluating polygon position against plane
 *	@tparam _TyVector is a vector type (specialization of Vector2 or Vector3)
 */
template <class _TyVector>
class CEvalPlanePos;

/**
 *	@brief small function object for evaluating polygon position against plane
 *		(specialization for 3D vertices)
 *	@tparam _TyScalar is a scalar type
 */
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER >= 1400
template <class _TyScalar>
class CEvalPlanePos<Vector3<_TyScalar> > {
#else // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
template <>
class CEvalPlanePos<Vector3f> { // MSVC6 doesn't support advanced template specialization
	typedef float _TyScalar;
#endif // !_MSC_VER || __MWERKS__ || _MSC_VER >= 1400
protected:
	size_t m_n_front_num; /**< @brief number of vertices in front of the plane */
	size_t m_n_back_num; /**< @brief number of vertices in back of the plane */
	const Plane3<_TyScalar> &m_r_t_plane; /**< @brief the classification plane */
	_TyScalar m_f_epsilon_ex; /**< @brief half of the plane thickness */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] r_t_plane is the classification plane
	 *	@param[in] f_epsilon_ex is half of the plane thickness
	 */
	inline CEvalPlanePos(const Plane3<_TyScalar> &r_t_plane, _TyScalar f_epsilon_ex)
		:m_r_t_plane(r_t_plane), m_n_front_num(0), m_n_back_num(0), m_f_epsilon_ex(f_epsilon_ex)
	{}

	/**
	 *	@brief inspects a polygon vertex
	 *	@param[in] r_v_vertex is a polygon vertex (order of vertices does not matter)
	 */
	inline void operator ()(const Vector3<_TyScalar> &r_v_vertex)
	{
		EPlanePos n_pos;
		if((n_pos = m_r_t_plane.n_Vector_Pos(r_v_vertex, m_f_epsilon_ex)) == plane_Front)
			++ m_n_front_num;
		else if(n_pos == plane_Back)
			++ m_n_back_num;
		// determine vertex side and increment counters accordingly
	}

	/**
	 *	@brief calculates position of the polygon relative to the plane
	 *	@return Returns the position of the given polygon relative to the given plane
	 *		(one of plane_Front).
	 *	@note This is untested - test it.
	 */
	inline operator EPlanePos() const
	{
		if(m_n_front_num) {
			if(!m_n_back_num)
				return plane_Front; // front propably occurres more frequently than split
			else
				return plane_Split;
		} else {
			if(m_n_back_num)
				return plane_Back; // back propably occurres more frequently than onplane
			else
				return plane_Onplane;
		}
	}
};

/**
 *	@brief a small class that finds a projection interval of polygon vertices onto an axis
 *	@tparam _TyVector is a vector type (specialization of Vector2 or Vector3)
 */
template <class _TyVector>
class CFindMinMaxProjection {
protected:
	typedef typename _TyVector::_TyScalar _TyScalar; /**< @brief a scalar type */
	_TyScalar m_f_min_projection; /**< @brief a minimal projection value */
	_TyScalar m_f_max_projection; /**< @brief a maximal projection value */
	_TyVector m_v_axis_dir; /**< @brief projection axis direction */
	_TyVector m_v_axis_org; /**< @brief point on the projection axis */

public:
	/**
	 *	@brief default constructor
	 *
	 *	@param[in] v_axis_dir is projection axis direction
	 *	@param[in] v_axis_org is point on the projection axis
	 *	@param[in] v_first_vertex is the first vertex to find projection for
	 */
	inline CFindMinMaxProjection(const _TyVector &v_axis_dir,
		const _TyVector &v_axis_org, const _TyVector &v_first_vertex)
		:m_f_min_projection(v_axis_dir.f_Dot(v_first_vertex - v_axis_org)),
		m_v_axis_dir(v_axis_dir), m_v_axis_org(v_axis_org)
	{
		m_f_max_projection = m_f_min_projection;
	}

	/**
	 *	@brief inspects a polygon vertex
	 *	@param[in] r_v_vertex is a polygon vertex (order of vertices does not matter)
	 */
	inline void operator ()(const Vector3<_TyScalar> &r_v_vertex)
	{
		_TyScalar f_projection = m_v_axis_dir.f_Dot(r_v_vertex - m_v_axis_org);
		if(m_f_min_projection > f_projection)
			m_f_min_projection = f_projection;
		if(m_f_max_projection < f_projection)
			m_f_max_projection = f_projection;
	}

	/**
	 *	@brief gets the minimal projection
	 *	@return Returns minimal projection, in multiples of projection axis vector.
	 */
	inline _TyScalar f_Min_Projection() const
	{
		return m_f_min_projection;
	}

	/**
	 *	@brief gets the maximal projection
	 *	@return Returns maximal projection, in multiples of projection axis vector.
	 */
	inline _TyScalar f_Max_Projection() const
	{
		return m_f_max_projection;
	}
};

/**
 *	@brief information about ray-triangle intersection
 */
struct TIntersectInfo_NoInfo {
	/**
	 *	@brief default constructor (has no effect)
	 */
	inline TIntersectInfo_NoInfo()
	{}

	/**
	 *	@brief constructor
	 *
	 *	@param[in] v_edge2 is the second edge of the triangle
	 *	@param[in] v_q is the q vector
	 *	@param[in] f_inv_det is the inverse determinant
	 *	@param[in] f_u is the first barycentric coordinate of the intersection
	 *	@param[in] f_v is the second barycentric coordinate of the intersection
	 *	@param[in] n_triangle is index of triangle where the intersection occured
	 */
	inline TIntersectInfo_NoInfo(const Vector3f &UNUSED(v_edge2), const Vector3f &UNUSED(v_q),
		float UNUSED(f_inv_det), float UNUSED(f_u), float UNUSED(f_v), size_t UNUSED(n_triangle))
	{}
};

/**
 *	@brief information about ray-triangle intersection
 */
struct TIntersectInfo_Distance {
	float f_t; /**< @brief distance from the intersection */

	/**
	 *	@brief default constructor (has no effect)
	 */
	inline TIntersectInfo_Distance()
	{}

	/**
	 *	@brief constructor
	 *
	 *	@param[in] v_edge2 is the second edge of the triangle
	 *	@param[in] v_q is the q vector
	 *	@param[in] f_inv_det is the inverse determinant
	 *	@param[in] f_u is the first barycentric coordinate of the intersection
	 *	@param[in] f_v is the second barycentric coordinate of the intersection
	 *	@param[in] n_triangle is index of triangle where the intersection occured
	 */
	inline TIntersectInfo_Distance(const Vector3f v_edge2, const Vector3f v_q,
		float f_inv_det, float UNUSED(f_u), float UNUSED(f_v), size_t UNUSED(n_triangle))
		:f_t(v_edge2.f_Dot(v_q) * f_inv_det)
	{}
};

/**
 *	@brief information about ray-triangle intersection
 */
struct TIntersectInfo_Full {
	float f_u; /**< @brief the first barycentric coordinate of the intersection */
	float f_v; /**< @brief the second barycentric coordinate of the intersection */
	size_t n_triangle; /**< @brief index of triangle where the intersection occured */
	float f_t; /**< @brief distance from the intersection */

	/**
	 *	@brief default constructor (has no effect)
	 */
	inline TIntersectInfo_Full()
	{}

	/**
	 *	@brief constructor
	 *
	 *	@param[in] v_edge2 is the second edge of the triangle
	 *	@param[in] v_q is the q vector
	 *	@param[in] f_inv_det is the inverse determinant
	 *	@param[in] _f_u is the first barycentric coordinate of the intersection
	 *	@param[in] _f_v is the second barycentric coordinate of the intersection
	 *	@param[in] _n_triangle is index of triangle where the intersection occured
	 */
	inline TIntersectInfo_Full(const Vector3f v_edge2, const Vector3f v_q,
		float f_inv_det, float _f_u, float _f_v, size_t _n_triangle)
		:f_u(_f_u), f_v(_f_v), n_triangle(_n_triangle),
		f_t(v_edge2.f_Dot(v_q) * f_inv_det)
	{}
};

/**
 *	@brief an empty type for use as void surface flags or material id
 */
struct CPolygon_EmptyType {
	/**
	 *	@brief swap two empty types (avoids unused variable warning in MSVC60)
	 *	@param[in] other is the other empty type to swap with (unused)
	 */
	inline void swap(CPolygon_EmptyType UNUSED(other))
	{}
};

/**
 *	@brief vertex traits class
 *	@tparam TVertStruct is type of vertex
 */
template <class TVertStruct>
class CVertexTraits {
public:
	typedef typename TVertStruct::TIntermediate TIntermediate; /**< @brief intermediate type for storing vertex position */
};

/**
 *	@brief vertex traits class (specialization for Vector2f)
 */
template <>
class CVertexTraits<Vector2f> {
public:
	typedef Vector2f TIntermediate; /**< @brief intermediate type for storing vertex position */
};

/**
 *	@brief vertex traits class (specialization for Vector3f)
 */
template <>
class CVertexTraits<Vector3f> {
public:
	typedef Vector3f TIntermediate; /**< @brief intermediate type for storing vertex position */
};

/**
 *	@brief vertex traits class (specialization for Vector4f)
 */
template <>
class CVertexTraits<Vector4f> {
public:
	typedef Vector4f TIntermediate; /**< @brief intermediate type for storing vertex position */
};

// todo - non-float vertices

} // ~poly2

/**
 *	@brief swap two empty types (avoids unused variable warning in MSVC60)
 *	@param[in] a is the first empty object to swap (unused)
 *	@param[in] b is the second empty object to swap (unused)
 */
inline void swap(poly2::CPolygon_EmptyType UNUSED(a),
	poly2::CPolygon_EmptyType UNUSED(b))
{}

namespace std {

/**
 *	@brief swap two empty types (avoids unused variable warning in MSVC60)
 *	@param[in] a is the first empty object to swap (unused)
 *	@param[in] b is the second empty object to swap (unused)
 */
inline void swap(poly2::CPolygon_EmptyType UNUSED(a),
	poly2::CPolygon_EmptyType UNUSED(b))
{}

} // ~std

/**
 *	@brief a simple (convex) polygon template
 *
 *	@tparam TVertStruct is a vertex type; needs to be convertible to Vector3f
 *	@tparam TSurfaceFlags is surface flags type (typically an integer or a void type)
 *	@tparam TMaterialId is material id type (typically an integer or a void type)
 */
template <class TVertStruct, class TSurfaceFlags = poly2::CPolygon_EmptyType,
	class TMaterialId = poly2::CPolygon_EmptyType>
class CPolygon2 {
public:
	typedef TVertStruct _TyVertex; /**< @brief vertex type of this polygon */
	typedef Vector3f _TyVector; /**< @brief position type for a given vertex type */
	typedef Plane3f _TyPlane; /**< @brief plane type for a given vertex type */
	typedef CPolygon2<TVertStruct, TSurfaceFlags, TMaterialId> _TyPolygon; /**< @brioef type of this template instance */

protected:
	std::vector<TVertStruct> m_vertex_list; /**< @brief list of polygon vertices */
	_TyPlane m_t_normal; /**< @brief polygon normal */

	// todo make some traits, support 2d polygons as well?

	TSurfaceFlags m_n_surface_flags; /**< @brief surface flags */
	TMaterialId m_n_material_id; /**< @brief material id */

public:
	/**
	 *	@brief default constructor (has no effect)
	 */
	inline CPolygon2()
	{}

	/**
	 *	@brief copy-constructor
	 *	@note This function throws std::bad_alloc.
	 */
	CPolygon2(const _TyPolygon &r_other) // throw std::bad_alloc)
		:m_t_normal(r_other.m_t_normal), m_n_surface_flags(r_other.m_n_surface_flags),
		m_n_material_id(r_other.m_n_material_id)
	{
		m_vertex_list.insert(m_vertex_list.begin(),
			r_other.m_vertex_list.begin(), r_other.m_vertex_list.end());
	}

	/**
	 *	@brief copy-operator
	 *	@note This function throws std::bad_alloc.
	 */
	_TyPolygon &operator =(const _TyPolygon &r_other) // throw(std::bad_alloc)
	{
		m_t_normal = r_other.m_t_normal;
		// copy normal

		m_n_surface_flags = r_other.m_n_surface_flags;
		m_n_material_id = r_other.m_n_material_id;
		// copy flags

		m_vertex_list.clear(); // avoid copying data
		m_vertex_list.insert(m_vertex_list.begin(),
			r_other.m_vertex_list.begin(), r_other.m_vertex_list.end());
		// copy vertex list (throws)

		return *this;
	}

	/**
	 *	@brief swaps two polygons
	 *	@param[in,out] r_t_other is the other polygon to swap with
	 */
	void Swap(_TyPolygon &r_t_other)
	{
		std::swap(m_t_normal, r_t_other.m_t_normal);
		std::swap(m_vertex_list, r_t_other.m_vertex_list);
		std::swap(m_n_surface_flags, r_t_other.m_n_surface_flags);
		std::swap(m_n_material_id, r_t_other.m_n_material_id);
	}

	/**
	 *	@brief deletes all polygon vertices
	 */
	inline void Delete()
	{
		m_vertex_list.clear();
	}

	/**
	 *	@brief inserts a single vertex at a specified location
	 *
	 *	@param[in] n_insert_before is a (zero-based) index of the vertex to insert the new vertices before
	 *	@param[in] r_t_vertex is a new vertex to be inserted
	 *
	 *	@note This function throws std::bad_alloc.
	 */
	inline void Insert_Vertex(size_t n_insert_before, const TVertStruct &r_t_vertex) // throw(std::bad_alloc)
	{
		m_vertex_list.insert(m_vertex_list.begin() + n_insert_before, r_t_vertex);
	}

	/**
	 *	@brief inserts a single vertex at the end of the vertex list
	 *	@param[in] r_t_vertex is a new vertex to be inserted
	 *	@note This function throws std::bad_alloc.
	 */
	inline void Insert_Vertex(const TVertStruct &r_t_vertex) // throw(std::bad_alloc)
	{
		m_vertex_list.insert(m_vertex_list.end(), r_t_vertex);
	}

	/**
	 *	@brief inserts several vertices at a specified location
	 *
	 *	@param[in] n_insert_before is a (zero-based) index of the vertex to insert the new vertices before
	 *	@param[in] p_vertex_it is const iterator to the first vertex to be inserted
	 *	@param[in] p_last_it is const iterator pointing to one past the last vertex to be inserted
	 *
	 *	@note This function throws std::bad_alloc.
	 */
	template <class CVertexIterator>
	inline void Insert_Vertex(size_t n_insert_before,
		CVertexIterator p_vertex_it, CVertexIterator p_last_it) // throw(std::bad_alloc)
	{
		m_vertex_list.insert(m_vertex_list.begin() + n_insert_before,
			p_vertex_it, p_last_it);
	}

	/**
	 *	@brief inserts several vertices at the end of the vertex list
	 *
	 *	@param[in] p_vertex_it is const iterator to the first vertex to be inserted
	 *	@param[in] p_last_it is const iterator pointing to one past the last vertex to be inserted
	 *
	 *	@note This function throws std::bad_alloc.
	 */
	template <class CVertexIterator>
	inline void Insert_Vertex(CVertexIterator p_vertex_it,
		CVertexIterator p_last_it) // throw(std::bad_alloc)
	{
		m_vertex_list.insert(m_vertex_list.end(), p_vertex_it, p_last_it);
	}

	/**
	 *	@brief deletes selected range of vertices
	 *
	 *	@param[in] n_index is (zero-based) index of the first vertex to be deleted
	 *	@param[in] n_vertex_num is number of vertices to be deleted
	 */
	inline void Delete_Vertices(size_t n_index, size_t n_vertex_num)
	{
		m_vertex_list.erase(m_vertex_list.begin() + n_index,
			m_vertex_list.begin() + (n_index + n_vertex_num));
	}

	/**
	 *	@brief gets number of vertices
	 *	@return Returns number of polygon vertices.
	 */
	inline size_t n_Vertex_Num() const
	{
		return m_vertex_list.size();
	}

	/**
	 *	@brief determines if number of vertices is nonzero
	 *	@return Returns true if number of polygon vertices is zero, otherwise returns false.
	 */
	inline bool b_Empty() const
	{
		_ASSERTE(m_vertex_list.empty() || m_vertex_list.size() >= 3); // make sure it is a valid polygon (at least triangle) if not empty
		return m_vertex_list.empty();
	}

	/**
	 *	@brief vertex access function
	 *	@param[in] n_index is zero-based vertex index
	 *	@return Returns reference to the selected vertex.
	 */
	inline TVertStruct &t_Vertex(size_t n_index)
	{
		return m_vertex_list[n_index];
	}

	/**
	 *	@brief vertex access function
	 *	@param[in] n_index is zero-based vertex index
	 *	@return Returns const reference to the selected vertex.
	 */
	inline const TVertStruct &t_Vertex(size_t n_index) const
	{
		return m_vertex_list[n_index];
	}

	/**
	 *	@brief gets surface flags
	 *	@return Returns polygon surface flags.
	 */
	inline TSurfaceFlags n_SurfaceFlags() const
	{
		return m_n_surface_flags;
	}

	/**
	 *	@brief gets material id
	 *	@return Returns polygon material id.
	 */
	inline TMaterialId n_MaterialId() const
	{
		return m_n_material_id;
	}

	/**
	 *	@brief sets surface flags
	 *	@param[in] n_surface_flags is the new value of surface flags
	 */
	inline void Set_SurfaceFlags(TSurfaceFlags n_surface_flags)
	{
		m_n_surface_flags = n_surface_flags;
	}

	/**
	 *	@brief sets material id
	 *	@param[in] n_material_id is the new value of material id
	 */
	inline void Set_MaterialId(TMaterialId n_material_id)
	{
		m_n_material_id = n_material_id;
	}

	/**
	 *	@brief calculates polygon normal plane
	 *	@return Return true on success, false if there were no three vertices to form a base.
	 *	@note This has to be called explicitly, it is never called by CPolygon2 itself.
	 */
	bool Calc_Normal(float f_epsilon_ex = f_epsilon)
	{
		// this only works for 3d anyway, probably need to template away

		if(m_vertex_list.size() < 3) {
			m_t_normal = Plane3f(Vector3f(0, 0, 0), 0);
			return false;
		}
		// impossible with so little vertices

#if 1
		Vector3f v_normal = ((Vector3f)m_vertex_list.back()).v_Cross((Vector3f)m_vertex_list.front());
		for(size_t i = 1, n = m_vertex_list.size(); i < n; ++ i)
			v_normal += ((Vector3f)m_vertex_list[i - 1]).v_Cross((Vector3f)m_vertex_list[i]);
		float f_length = v_normal.f_Length();
		if(f_length >= f_epsilon_ex) {
			v_normal /= -f_length; // note the minus; the normal calculated by the algorithm below is negative
			m_t_normal = _TyPlane((Vector3f)m_vertex_list.front(), v_normal, plane::from_position_normal);
			// create the plane
			return true;
		}
		// use Newell's algorithm, this is probably more robust (although still ...)
#else // 1
		float f_short_normal_length = 0;
		Vector3f v_short_normal;

		const float f_epsilon_ex2 = f_epsilon_ex * f_epsilon_ex;
		for(size_t i = 1, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_u = (Vector3f)m_vertex_list[i - 1] - (Vector3f)m_vertex_list[i];
			float f_ulen2;
			if((f_ulen2 = v_u.f_Length2()) < f_epsilon_ex2)
				continue;
			// get a non-zero vector

			size_t b = i + 1, e = n; // try all vertices after i
			for(size_t n_pass = 0; n_pass < 2; ++ n_pass) { // unroll?
				for(size_t j = b; j < e; ++ j) {
					_ASSERTE(i != j && j != i - 1);
					// make sure we're not wasting time on degenerate cases

					Vector3f v_v = (Vector3f)m_vertex_list[j] - (Vector3f)m_vertex_list[i];
					float f_vlen2;
					if((f_vlen2 = v_v.f_Length2()) < f_epsilon_ex2)
						continue;
					// get another different non-zero vector

#if 1
					Vector3f v_normal = v_u.v_Cross(v_v);
					// for unnormalized vectors, the magnitude is ||u|| * ||v|| * sin(theta), if both vectors
					// are epsilon length, it will be at most f_epsilon_ex2. normalization is imperative.

					float f_normalized_length2 = v_normal.f_Length2();
					float f_inv_length2 = 1 / f_normalized_length2;
					f_normalized_length2 /= (f_ulen2 * f_vlen2);
					// could divide the whole magnitude by sqrt(f_vlen2 * f_ulen2)
#else // 1
					Vector3f v_normal = (v_u / float(sqrt(f_vlen2))).v_Cross(v_v / float(sqrt(f_ulen2)));
					// normalizing the vectors themselves is more expensive, maybe slightly more numerically precise

					float f_normalized_length2 = v_normal.f_Length2();
					float f_inv_length2 = 1 / f_normalized_length2;
#endif // 1
					if(f_normalized_length2 < 1e-3f/*f_epsilon_ex2*/) { // amounts to sin(theta)^2, for epsilon = 0.001 it is about 0.06 degrees, for epsilon sqrt(0.001) it is about 2 debrees. thresholding too low here will lead to accepting shallow edges as normals, and increase noise in the normals. the normals should be ranked and the longest (the best defined) should be returned. this greedy approach is not very robust.
						if(f_short_normal_length < f_normalized_length2) {
							f_short_normal_length = f_normalized_length2;
							v_short_normal = v_normal;
						}
						continue;
					}
					v_normal *= sqrt(f_inv_length2); // normalise using the inverse
					_ASSERTE(fabs(v_normal.f_Length() - 1) < 1e-3f);
					// colinear vectors would yield zero normal

					m_t_normal = _TyPlane((Vector3f)m_vertex_list[i], v_normal, plane::from_position_normal);
					// create the plane

					return true;
					// success: found a base that is not degenerate and can be orthogonalized
				}

				b = 0;
				e = max(size_t(2), i) - 2;
				// try vertices before i (there may be none if i == 1)
			}
		}
		// try to come up with a (reasonably) strong normal

		// note that this loop may still give poor results: in case the polygon is a finely subdivided
		// circle, the edges may be very short (it will take one short edge i-1, i and one long i, j that will be sufficiently ), but overall, the plane is defined very well
		// it would probably be better to do some SVD (with all the vertices) in these cases.
		// let's hope we will mostly have nice small polygons

#if 1
		if(f_short_normal_length > 1e-6f) {
			v_short_normal.Normalize();
			m_t_normal = _TyPlane((Vector3f)m_vertex_list[0], v_short_normal, plane::from_position_normal);
			// create the plane

			return true;
			// success: found a base that is not degenerate and can be orthogonalized
		}
		// resort to using a shorter plane 
#endif // 1
#endif // 1

		m_t_normal = Plane3f(Vector3f(0, 0, 0), 0);

		return false;
	}

	/**
	 *	@brief normal access function
	 *	@return Returns polygon normal.
	 */
	inline const _TyPlane &t_Normal() const
	{
		return m_t_normal;
	}

	/**
	 *	@brief for-each cycle over polygon's vertices
	 *	@tparam CFunctionObject is a function object type
	 *	@param[in] t_object is an instance of the function object
	 *	@return Returns the function object, after passing all the vertices to it.
	 */
	template <class CFunctionObject>
	inline CFunctionObject for_each_vertex(CFunctionObject t_object)
	{
		return std::for_each(m_vertex_list.begin(), m_vertex_list.end(), t_object);
	}

	/**
	 *	@brief for-each cycle over polygon's (const) vertices
	 *	@tparam CFunctionObject is a function object type
	 *	@param[in] t_object is an instance of the function object
	 *	@return Returns the function object, after passing all the vertices to it.
	 */
	template <class CFunctionObject>
	inline CFunctionObject for_each_vertex(CFunctionObject t_object) const
	{
		return std::for_each(m_vertex_list.begin(), m_vertex_list.end(), t_object);
	}

	/**
	 *	@brief swaps vertex order, flips normal
	 */
	void Swap_VertexOrder()
	{
		m_t_normal.v_normal *= -1;
		m_t_normal.f_dist *= -1;

		for(size_t i = 0, n = m_vertex_list.size() / 2,
		   j = m_vertex_list.size() - 1; i < n; ++ i, -- j)
			std::swap(m_vertex_list[i], m_vertex_list[j]);
		// swap vertices
	}

	/**
	 *	@brief gets position of the polygon center
	 *	@return Returns the position of the center of this polygon,
	 *		if it has no vertices, returns the O vector.
	 */
	const _TyVector v_Center() const
	{
		if(m_vertex_list.empty()) {
			_TyVector v;
			v.SetZero();
			return v;
		}

		_TyVector v_center;
		v_center.SetZero();
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i)
			v_center += (_TyVector)m_vertex_list[i];

		return v_center * (1.0f / m_vertex_list.size());
	}

	/**
	 *	@brief gets coordinates of the polygon center
	 *	@return Returns the position of the center of this polygon,
	 *		if it has no vertices, returns the O vertex.
	 */
	const typename poly2::CVertexTraits<TVertStruct>::TIntermediate t_Center() const // works with new MSVC, will probably fail on MSVC60 (but MSVC 60 doeesn't perform much checking, can simply replace traits by TVertStruct::TIntermediate and don't call this function with Vector3f)
	{
		_TyVector v_center;
		v_center.SetZero();
		typename poly2::CVertexTraits<TVertStruct>::TIntermediate t_center(v_center);
		if(m_vertex_list.empty())
			return t_center;

		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i)
			t_center += m_vertex_list[i];

		return t_center * (1.0f / m_vertex_list.size());
	}

	/**
	 *	@brief polygon convexness test
	 *	@param[in] f_epsilon_ex is maximal tolerated concave edge height
	 *		(negative values will reject convex polygons as well)
	 *	@return Returns true in case the polygon is convex (or empty), otherwise returns false.
	 *	@note This function requires the polygon normal to be up to date.
	 */
	bool b_IsConvex(float f_epsilon_ex = f_epsilon) const
	{
		size_t n;
		if((n = m_vertex_list.size()) < 3)
			return true;
		// polygon with so little vertices can't be concave

		return std::find_if(m_vertex_list.begin(), m_vertex_list.end(),
			poly2::CFindConcave<_TyVector>(m_vertex_list[n - 2], m_vertex_list[n - 1],
			m_t_normal.v_normal, f_epsilon_ex)) == m_vertex_list.end();
		// in case concave edge is not found, it's convex polygon
	}

	/**
	 *	@brief gets polygon plane deviation
	 *	@return Returns maximal absolute distance of polygon points from polygon plane.
	 *	@note This function requires the polygon normal to be up to date.
	 */
	inline float f_MaxPlaneDeviation() const
	{
		return std::for_each(m_vertex_list.begin(), m_vertex_list.end(),
			poly2::CFindMaxPlaneDeviation<_TyVector>(m_t_normal));
	}

	/**
	 *	@brief determines whether the polygon is tiny
	 *	@param[in] f_epsilon_ex is tiny edge length threshold
	 *	@return Returns true in case polygon doesn't contain at
	 *		least three edges longer than epsilon, otherwise returns false.
	 */
	inline bool b_IsTiny(float f_epsilon_ex = f_epsilon) const
	{
		if(m_vertex_list.size() < 3)
			return true;
		return std::find_if(m_vertex_list.begin(), m_vertex_list.end(),
			poly2::CFindLongEdges<_TyVector>(m_vertex_list.back(), f_epsilon_ex)) == m_vertex_list.end();
	}

	/**
	 *	@brief calculates position of the polygon relative to a plane
	 *
	 *	@param[in] r_t_plane is the plane
	 *	@param[in] f_epsilon_ex is half of the plane thickness
	 *
	 *	@return Returns position of the polygon relative to the given plane
	 *		(one of plane_Front, plane_Back, plane_Split or plane_Onplane).
	 *
	 *	@note Polygons with no vertices are always classified as plane_Onplane.
	 */
	inline EPlanePos n_Plane_Pos(const _TyPlane &r_t_plane, float f_epsilon_ex = f_epsilon) const
	{
		return std::for_each(m_vertex_list.begin(),
			m_vertex_list.end(), poly2::CEvalPlanePos<_TyVector>(r_t_plane, f_epsilon_ex));
	}

	/**
	 *	@brief tests if a point is inside the polygon
	 *
	 *	@param[in] v_point is the tested point
	 *	@param[in] f_epsilon_ex is maximal angular discrepancy to pass the test (in radians)
	 *
	 *	@return Returns true in case the given point is inside the polygon.
	 *
	 *	@note This function uses method of summing up point-to-edge vector angles.
	 */
	bool b_PointInside_SumAngles(Vector3f v_point, float f_epsilon_ex = f_epsilon) const
	{
		if(m_vertex_list.empty())
			return false;

		float f_angle_sum = 0;
		Vector3f v_prev = (Vector3f)m_vertex_list.back() - v_point;
		v_prev.Normalize();
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_cur = (Vector3f)m_vertex_list[i] - v_point;
			v_cur.Normalize();
			// take two vectors ...

			f_angle_sum += (float)acos(v_prev.f_Dot(v_cur));
			// sum up the angles they hold

			v_prev = v_cur;
		}
		// sum up the angles

		return fabs(f_angle_sum - 2 * f_pi) < f_epsilon_ex;
		// in case difference from 2pi is less than epsilon, the point is inside
	}

	/**
	 *	@brief tests if a point is inside the polygon (or above it)
	 *
	 *	@param[in] v_point is the tested point
	 *	@param[in] f_epsilon_ex is maximal distance of the point from the planes
	 *		(bigger values makes more points pass the test)
	 *
	 *	@return Returns true in case the given point is inside the polygon.
	 *
	 *	@note This function requires the polygon normal to be up to date.
	 *	@note This function uses method of edge planes 'fence'.
	 *	@note This function accepts all points that are inside the polygon after projection
	 *		on it (it does not check if the point actually lies in the polygon plane).
	 */
	bool b_PointInside_EdgePlanes(Vector3f v_point, float f_epsilon_ex = f_epsilon) const
	{
		Vector3f v_prev = (Vector3f)m_vertex_list.back();
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_cur = (Vector3f)m_vertex_list[i];
			Vector3f v_u = v_prev - v_cur;
			v_prev = v_cur;
			// get an edge vector ...

			Vector3f v_v = v_point - v_cur;
			// get vector to point

			if(v_u.v_Cross(m_t_normal.v_normal).f_Dot(v_v) > f_epsilon)
				return false;
			// see if the point lies behind the 'fence' plane
			// (in case it does, it's outside)
		}

		return true;
	}

	/**
	 *	@brief calculates an edge plane (the polygon is on its front side)
	 *
	 *	@param[in] n_edge is zero-based index of an edge
	 *	@param[in] f_offset is offset of the plane along its normal
	 *		(positive offset pushes the plane inside the polygon)
	 *
	 *	@note This function requires the polygon normal to be up to date.
	 *	@note In case the polygon is not convex, sidedness of the normal is undefined.
	 *
	 *	@todo This is untested - test it.
	 */
	Plane3f t_EdgePlane(size_t n_edge, float f_offset = 0) const
	{
		_ASSERTE(!m_vertex_list.empty());
		size_t n = m_vertex_list.size();
		Vector3f v_prev = (Vector3f)m_vertex_list[n_edge];
		Vector3f v_cur = (Vector3f)m_vertex_list[(n_edge + 1) % n];
		Vector3f v_edge_normal = (v_cur - v_prev).v_Cross(m_t_normal.v_normal);
		v_edge_normal.Normalize();
		Plane3f t_plane(v_prev, v_edge_normal, plane::from_position_normal);

		_ASSERTE(!b_IsConvex() || t_plane.n_Vector_Pos(v_Center()) == plane_Front);
		// make sure that the center is in front of the plane (must not be true for concave polygons)

		t_plane.f_dist -= f_offset;

		return t_plane;
	}

	/**
	 *	@brief tests if a point is inside the polygon
	 *
	 *	@param[in] v_point is the tested point
	 *	@param[in] f_epsilon_ex is maximal barycentric coordinate deviation
	 *		(bigger values makes more points pass the test)
	 *
	 *	@return Returns true in case the given point is inside the polygon.
	 *
	 *	@note This function uses barycentric coordinates and polygon to triangles decomposition.
	 *	@note In here, the epsilon is relative to the poslygon size, unlike in the other versions of the test.
	 */
	bool b_PointInside_Barycentric(Vector3f v_point, float f_epsilon_ex = 0) const
	{
		Vector3f v_vert0 = (Vector3f)m_vertex_list.front();
		Vector3f v_vert1 = (Vector3f)m_vertex_list[1];
		for(size_t i = 2, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_vert2 = (Vector3f)m_vertex_list[i];

			Vector3f v_edge1 = v_vert1 - v_vert0,
					 v_edge2 = v_vert2 - v_vert0;
			v_vert1 = v_vert2; // not used below
			float f_inv_length2_edge1 = 1.0f / v_edge1.f_Length2();
			Vector3f v_perp = v_edge2 - v_edge1 * (v_edge1.f_Dot(v_edge2) * f_inv_length2_edge1);
			Vector3f v_t = v_point - v_vert0;
			float f_u = v_t.f_Dot(v_perp) / v_perp.f_Dot(v_edge2);
			// calculate barycentric coordinates for v_point

			if(f_u < -f_epsilon_ex)
				continue;
			// first condition failed

			float f_v = (v_t - v_edge2 * f_u).f_Dot(v_edge1) * f_inv_length2_edge1;
			if(f_v >= -f_epsilon_ex && f_u + f_v <= 1 + f_epsilon_ex)
				return true;
			// in case the other conditions are satisfied as well, point lies inside
		}

		return false;
	}

	/**
	 *	@brief determines if two polygons overlap, using (3D) separating axes theorem
	 *
	 *	@tparam CPolygon2_TVertStruct2_ is a second polgon type
	 *
	 *	@param[in] r_poly2 is another polygon
	 *	@param[out] r_n_sep_axis is maximum separation axis (zero-based) index
	 *		(points to the second vertex of the edge, -1 means separated by normal plane)
	 *	@param[out] r_b_sep_axis_on_poly2 is separating axis location flag; if set,
	 *		it is on r_poly2, otherwise it is on this
	 *	@param[out] r_f_sep_distance is separating distance (if below epsilon, overlap occurs)
	 *	@param[in] f_epsilon_ex is epsilon for the collision
	 *		(maximal allowed projection overlap; greater values cause more overlaps)
	 *
	 *	@return Returns true if polygons overlap, otherwise returns false.
	 *
	 *	@note This function requires the polygon normal to be up to date.
	 */
	template <class CPolygon2_TVertStruct2_>
	bool b_Overlap_SepAxes(const CPolygon2_TVertStruct2_ &r_poly2, size_t &r_n_sep_axis,
		bool &r_b_sep_axis_on_poly2, float &r_f_sep_distance, float f_epsilon_ex = f_epsilon) const
	{
		if(m_vertex_list.empty() || r_poly2.b_Empty()) {
			r_n_sep_axis = size_t(-1);
			r_b_sep_axis_on_poly2 = false;
			r_f_sep_distance = 0; // ???
			return false; // nothing to overlap
		}

		size_t n_sep_axis;
		float f_max_distance;
		{
			Vector3f v_vert0 = (Vector3f)m_vertex_list.back();

			poly2::CFindMinMaxProjection<_TyVector> proj(m_t_normal.v_normal, v_vert0, m_vertex_list.front());
			proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj); // unnecessary if it's truly planar
			poly2::CFindMinMaxProjection<_TyVector> proj2(m_t_normal.v_normal, v_vert0, r_poly2.t_Vertex(0));
			proj2 = r_poly2.for_each_vertex(proj2);
			// test all the vertices agains the normal plane of this polygon (handles 3D cases)

			/*bool b_separation = proj2.f_Min_Projection() > proj.f_Max_Projection() ||
				proj.f_Min_Projection() > proj2.f_Max_Projection();*/ // unused
			// tell wheter axis is a separation axis (i.e. polygons do not collide on it)

			float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
				proj.f_Min_Projection() - proj2.f_Max_Projection());
			// find distance

			n_sep_axis = size_t(-1); // marks the normal plane
			f_max_distance = f_distance;
			// use this as the minimal separating distance

			for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
				Vector3f v_vert1 = (Vector3f)m_vertex_list[i];

				Vector3f v_org = v_vert0;
				Vector3f v_dir = m_t_normal.v_normal.v_Cross(v_vert1 - v_vert0);
				v_dir.Normalize();
				if(v_dir.f_Length2() < .5f)
					continue;
				// get a separating axis (needs polygon normal)

				poly2::CFindMinMaxProjection<_TyVector> proj(v_dir, v_org, m_vertex_list.front());
				proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
				poly2::CFindMinMaxProjection<_TyVector> proj2(v_dir, v_org, r_poly2.t_Vertex(0));
				proj2 = r_poly2.for_each_vertex(proj2);
				// test all the vertices

				/*bool b_separation = proj2.f_Min_Projection() > proj.f_Max_Projection() ||
					proj.f_Min_Projection() > proj2.f_Max_Projection();*/ // unused
				// tell wheter axis is a separation axis (i.e. polygons do not collide on it)

				float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
					proj.f_Min_Projection() - proj2.f_Max_Projection());
				// find distance

				if(f_max_distance < f_distance || n_sep_axis == size_t(-1)) {
					f_max_distance = f_distance;
					n_sep_axis = i;
				}
				// see if this axis has a bigger separating distance

				v_vert0 = v_vert1;
			}
			// test all the vertices of r_poly2 agains the edge planes of this polygon (handles 3D cases)
		}
		// calculates all separating axes for this polygon

		bool b_sep_axis_on_other = false;
		{
			Vector3f v_vert0 = (Vector3f)r_poly2.t_Vertex(r_poly2.n_Vertex_Num() - 1);

			poly2::CFindMinMaxProjection<_TyVector> proj(r_poly2.t_Normal().v_normal, v_vert0, m_vertex_list.front());
			proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
			poly2::CFindMinMaxProjection<_TyVector> proj2(r_poly2.t_Normal().v_normal, v_vert0, r_poly2.t_Vertex(0));
			proj2 = r_poly2.for_each_vertex(proj2); // unnecessary if it's truly planar
			// test all the vertices of r_poly2 agains the normal plane of this polygon (handles 3D cases)

			/*bool b_separation = proj2.f_Min_Projection() > proj.f_Max_Projection() ||
				proj.f_Min_Projection() > proj2.f_Max_Projection();*/ // unused
			// tell wheter axis is a separation axis (i.e. polygons do not collide on it)

			float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
				proj.f_Min_Projection() - proj2.f_Max_Projection());
			// find distance

			if(f_max_distance < f_distance || n_sep_axis == size_t(-1)) {
				f_max_distance = f_distance;
				n_sep_axis = size_t(-1); // marks the normal plane
				b_sep_axis_on_other = true;
			}
			// use this as the minimal separating distance

			for(size_t i = 0, n = r_poly2.n_Vertex_Num(); i < n; ++ i) {
				Vector3f v_vert1 = (Vector3f)r_poly2.t_Vertex(i);

				Vector3f v_org = v_vert0;
				Vector3f v_dir = r_poly2.t_Normal().v_normal.v_Cross(v_vert1 - v_vert0);
				v_dir.Normalize();
				if(v_dir.f_Length2() < .5f)
					continue;
				// get a separating axis (needs polygon normal)

				poly2::CFindMinMaxProjection<_TyVector> proj(v_dir, v_org, m_vertex_list.front());
				proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
				poly2::CFindMinMaxProjection<_TyVector> proj2(v_dir, v_org, r_poly2.t_Vertex(0));
				proj2 = r_poly2.for_each_vertex(proj2);
				// test all the vertices

				/*bool b_separation = proj2.f_Min_Projection() > proj.f_Max_Projection() ||
					proj.f_Min_Projection() > proj2.f_Max_Projection();*/ // unused
				// tell wheter axis is a separation axis (i.e. polygons do not collide on it)

				float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
					proj.f_Min_Projection() - proj2.f_Max_Projection());
				// find distance

				if(f_max_distance < f_distance || n_sep_axis == size_t(-1)) {
					f_max_distance = f_distance;
					n_sep_axis = i;
					b_sep_axis_on_other = true;
				}
				// see if this axis has a bigger separating distance

				v_vert0 = v_vert1;
			}
			// test all the vertices of r_poly2 agains the edge planes of this polygon (handles 3D cases)
		}
		// calculates all the separating axes for the other polygon

		r_n_sep_axis = n_sep_axis;
		r_b_sep_axis_on_poly2 = b_sep_axis_on_other;
		r_f_sep_distance = f_max_distance;
		// fill the outputs

		return f_max_distance < f_epsilon_ex; // in case maximum distance is below zero, there is overlap
	}

	/**
	 *	@brief determines if two polygons overlap, using (3D) separating axes theorem
	 *
	 *	@tparam CPolygon2_TVertStruct2_ is a second polgon type
	 *
	 *	@param[in] r_poly2 is another polygon
	 *	@param[in] f_epsilon_ex is epsilon for the collision
	 *		(maximal allowed projection overlap; greater values cause more overlaps)
	 *
	 *	@return Returns true if polygons overlap, otherwise returns false.
	 *
	 *	@note This function requires the polygon normal to be up to date.
	 */
	template <class CPolygon2_TVertStruct2_>
	bool b_Overlap_SepAxes(const CPolygon2_TVertStruct2_ &r_poly2, float f_epsilon_ex = f_epsilon) const
	{
		if(m_vertex_list.empty() || r_poly2.b_Empty())
			return false; // nothing to overlap

		{
			Vector3f v_vert0 = (Vector3f)m_vertex_list.back();

			poly2::CFindMinMaxProjection<_TyVector> proj(m_t_normal.v_normal, v_vert0, m_vertex_list.front());
			proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj); // unnecessary if it's truly planar
			poly2::CFindMinMaxProjection<_TyVector> proj2(m_t_normal.v_normal, v_vert0, r_poly2.t_Vertex(0));
			//proj2 = std::for_each(r_poly2.m_vertex_list.begin() + 1, r_poly2.m_vertex_list.end(), proj2); // unnecessary if it's truly planar
			proj2 = r_poly2.for_each_vertex(proj2); // loops once more than necessary - todo - overload, add indices
			// test all the vertices agains the normal plane of this polygon (handles 3D cases)

			float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
				proj.f_Min_Projection() - proj2.f_Max_Projection());
			// find distance

			if(f_distance >= f_epsilon_ex)
				return false;
			// there is a separating axis, no collision

			for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
				Vector3f v_vert1 = (Vector3f)m_vertex_list[i];

				Vector3f v_org = v_vert0;
				Vector3f v_dir = m_t_normal.v_normal.v_Cross(v_vert1 - v_vert0);
				v_dir.Normalize();
				if(v_dir.f_Length2() < .5f)
					continue;
				// get a separating axis (needs polygon normal)

				poly2::CFindMinMaxProjection<_TyVector> proj(v_dir, v_org, m_vertex_list.front());
				proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
				poly2::CFindMinMaxProjection<_TyVector> proj2(v_dir, v_org, r_poly2.t_Vertex(0));
				//proj2 = std::for_each(r_poly2.m_vertex_list.begin() + 1, r_poly2.m_vertex_list.end(), proj2);
				proj2 = r_poly2.for_each_vertex(proj2);
				// test all the vertices

				float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
					proj.f_Min_Projection() - proj2.f_Max_Projection());
				// find distance

				if(f_distance >= f_epsilon_ex)
					return false;
				// there is a separating axis, no collision

				v_vert0 = v_vert1;
			}
			// test all the vertices of r_poly2 agains the edge planes of this polygon (handles 3D cases)
		}
		// calculates all separating axes for this polygon

		{
			Vector3f v_vert0 = (Vector3f)r_poly2.t_Vertex(r_poly2.n_Vertex_Num() - 1);

			poly2::CFindMinMaxProjection<_TyVector> proj(r_poly2.t_Normal().v_normal, v_vert0, m_vertex_list.front());
			proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
			poly2::CFindMinMaxProjection<_TyVector> proj2(r_poly2.t_Normal().v_normal, v_vert0, r_poly2.t_Vertex(0));
			//proj2 = std::for_each(r_poly2.m_vertex_list.begin() + 1, r_poly2.m_vertex_list.end(), proj2); // unnecessary if it's truly planar
			proj2 = r_poly2.for_each_vertex(proj2); // unnecessary if it's truly planar
			// test all the vertices of r_poly2 agains the normal plane of this polygon (handles 3D cases)

			float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
				proj.f_Min_Projection() - proj2.f_Max_Projection());
			// find distance

			if(f_distance >= f_epsilon_ex)
				return false;
			// there is a separating axis, no collision

			for(size_t i = 0, n = r_poly2.n_Vertex_Num(); i < n; ++ i) {
				Vector3f v_vert1 = (Vector3f)r_poly2.t_Vertex(i);

				Vector3f v_org = v_vert0;
				Vector3f v_dir = r_poly2.t_Normal().v_normal.v_Cross(v_vert1 - v_vert0);
				v_dir.Normalize();
				if(v_dir.f_Length2() < .5f)
					continue;
				// get a separating axis (needs polygon normal)

				poly2::CFindMinMaxProjection<_TyVector> proj(v_dir, v_org, m_vertex_list.front());
				proj = std::for_each(m_vertex_list.begin() + 1, m_vertex_list.end(), proj);
				poly2::CFindMinMaxProjection<_TyVector> proj2(v_dir, v_org, r_poly2.t_Vertex(0));
				//proj2 = std::for_each(r_poly2.m_vertex_list.begin() + 1, r_poly2.m_vertex_list.end(), proj2);
				proj2 = r_poly2.for_each_vertex(proj2);
				// test all the vertices

				float f_distance = max(proj2.f_Min_Projection() - proj.f_Max_Projection(),
					proj.f_Min_Projection() - proj2.f_Max_Projection());
				// find distance

				if(f_distance >= f_epsilon_ex)
					return false;
				// there is a separating axis, no collision

				v_vert0 = v_vert1;
			}
			// test all the vertices of r_poly2 agains the edge planes of this polygon (handles 3D cases)
		}
		// calculates all the separating axes for the other polygon

		return true; // no separating axis found - have an overlap
	}

	/**
	 *	@brief performs a polygon-ray hit test
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_NoInfo t_no_info;
		return b_RayHit<false, false>(r_v_org, r_v_dir, t_no_info, f_epsilon_ex);
	}

	/**
	 *	@brief performs a polygon-ray hit test, culls front-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullFront(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_NoInfo t_no_info;
		return b_RayHit<true, false>(r_v_org, r_v_dir, t_no_info, f_epsilon_ex);
	}

	/**
	 *	@brief performs a polygon-ray hit test, culls back-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullBack(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_NoInfo t_no_info;
		return b_RayHit<false, true>(r_v_org, r_v_dir, t_no_info, f_epsilon_ex);
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit distance
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Distance t_dist;
		bool b_result = b_RayHit<false, false>(r_v_org, r_v_dir, t_dist, f_epsilon_ex);
		r_f_t = t_dist.f_t;
		return b_result;
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit distance, culls front-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullFront(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Distance t_dist;
		bool b_result = b_RayHit<true, false>(r_v_org, r_v_dir, t_dist, f_epsilon_ex);
		r_f_t = t_dist.f_t;
		return b_result;
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit distance, culls back-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullBack(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Distance t_dist;
		bool b_result = b_RayHit<poly2::TIntersectInfo_NoInfo,
			false, true>(r_v_org, r_v_dir, t_dist, f_epsilon_ex);
		r_f_t = t_dist.f_t;
		return b_result;
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit coordinates
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[out] r_f_u is the first barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_f_v is the second barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_n_triangle is zero-based index of the vertex forming the triangle with vertices
	 *		0 and 1 where the intersection occured (to be able to use the barycentric coordinates)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float &r_f_u, float &r_f_v, size_t &r_n_triangle,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Full t_iinfo;
		bool b_result = b_RayHit<poly2::TIntersectInfo_NoInfo,
			false, false>(r_v_org, r_v_dir, t_iinfo, f_epsilon_ex);
		r_f_t = t_iinfo.f_t;
		r_f_u = t_iinfo.f_u;
		r_f_v = t_iinfo.f_v;
		r_n_triangle = t_iinfo.n_triangle;
		return b_result;
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit coordinates,
	 *		culls front-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[out] r_f_u is the first barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_f_v is the second barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_n_triangle is zero-based index of the vertex forming the triangle with vertices
	 *		0 and 1 where the intersection occured (to be able to use the barycentric coordinates)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullFront(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float &r_f_u, float &r_f_v, size_t &r_n_triangle,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Full t_iinfo;
		bool b_result = b_RayHit<poly2::TIntersectInfo_NoInfo,
			true, false>(r_v_org, r_v_dir, t_iinfo, f_epsilon_ex);
		r_f_t = t_iinfo.f_t;
		r_f_u = t_iinfo.f_u;
		r_f_v = t_iinfo.f_v;
		r_n_triangle = t_iinfo.n_triangle;
		return b_result;
	}

	/**
	 *	@brief performs a polygon-ray hit test and calculates hit coordinates,
	 *		culls back-facing rays
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[out] r_f_t is hit distance, in multiples of r_v_dir (only valid if hit occured)
	 *	@param[out] r_f_u is the first barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_f_v is the second barycentric coordinate of the hit (only valid if hit occured)
	 *	@param[out] r_n_triangle is zero-based index of the vertex forming the triangle with vertices
	 *		0 and 1 where the intersection occured (to be able to use the barycentric coordinates)
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	inline bool b_RayHit_CullBack(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		float &r_f_t, float &r_f_u, float &r_f_v, size_t &r_n_triangle,
		float f_epsilon_ex = 1e-5f) const
	{
		poly2::TIntersectInfo_Full t_iinfo;
		bool b_result = b_RayHit<poly2::TIntersectInfo_NoInfo,
			false, true>(r_v_org, r_v_dir, t_iinfo, f_epsilon_ex);
		r_f_t = t_iinfo.f_t;
		r_f_u = t_iinfo.f_u;
		r_f_v = t_iinfo.f_v;
		r_n_triangle = t_iinfo.n_triangle;
		return b_result;
	}

	/**
	 *	@brief gets nearest point inside the polygon
	 *	@param[in] r_v_vector is a reference point
	 *	@return Returns point that is nearest to the given position and is inside the polygon.
	 */
	Vector3f v_NearestPoint_Barycentric(const Vector3f &r_v_vector) const
	{
		Vector3f v_projection = r_v_vector - m_t_normal.v_normal *
			(m_t_normal.v_normal.f_Dot(r_v_vector) + m_t_normal.f_dist);
		if(b_PointInside_Barycentric(v_projection, 0))
			return v_projection;
		// calculate projection, then in case it's inside, return it

		/*
		if(b_PointInside_EdgePlanes(r_v_vector)) {
			return r_v_vector - m_t_normal.v_normal *
				(m_t_normal.v_normal.f_Dot(r_v_vector) + m_t_normal.f_dist);
		}
		// in case it's inside, calculate projection and return it
		*/ // possibly faster version using edge planes 'fence'

		Vector3f v_closest;
		float f_min_dist2 = 1e10f; // actually any
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_u = (Vector3f)m_vertex_list[(i + 1) % n] - (Vector3f)m_vertex_list[i];
			// get an edge vector ...

			Vector3f v_v = r_v_vector - (Vector3f)m_vertex_list[i];
			// get vector from edge origin to point

			float f_t = v_u.f_Dot(v_v) / v_u.f_Length2();
			//float f_t = v_u.f_Dot(v_v) / v_u.f_Dot(v_u); // equivalent to the previous one
			// get projection distance along the edge

			Vector3f v_projection;
			if(f_t <= 0)
				v_projection = m_vertex_list[i];
			else if(f_t >= 1)
				v_projection = m_vertex_list[(i + 1) % n];
			else
				v_projection = (Vector3f)m_vertex_list[i] + v_u * f_t;
			// get projection position

			float f_dist2 = (r_v_vector - v_projection).f_Length2();
			if(f_dist2 < f_min_dist2 || !i) {
				f_min_dist2 = f_dist2;
				v_closest = v_projection;
			}
			// find the closest projection
		}

		return v_closest;
	}

	/**
	 *	@brief cuts the polygon by a plane so only the half lying on the specified side remains
	 *
	 *	@param[in] r_t_plane is the splitting plane
	 *	@param[in] n_desired_half is the desired half of the polygon to remain after
	 *		cutting (plane_Front and plane_Back are the only allowed values)
	 *	@param[in] f_epsilon_ex is half of the plane thickness (edges shorter
	 *		than that will not be split)
	 *
	 *	@note This function throws std::bad_alloc.
	 *	@note This function should not be used to cut a polygon in two,
	 *		it will not be surface-preserving. Use Split() instead.
	 */
	void Cut(const _TyPlane &r_t_plane, EPlanePos n_desired_half = plane_Front,
		float f_epsilon_ex = f_epsilon) // throw(std::bad_alloc)
	{
		_ASSERTE(n_desired_half == plane_Front || n_desired_half == plane_Back);

		if(m_vertex_list.empty())
			return;
		// nothing to do here

		EPlanePos n_other_half = (n_desired_half == plane_Front)? plane_Back : plane_Front;
		// the other ('wrong one') position

		TVertStruct t_prev_vertex = m_vertex_list.back();
		bool b_prev_good_side = r_t_plane.n_Vector_Pos(t_prev_vertex, f_epsilon_ex) != n_other_half;
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			TVertStruct t_vertex = m_vertex_list[i];
			bool b_good_side = r_t_plane.n_Vector_Pos(t_vertex, f_epsilon_ex) != n_other_half;
			// determine plane sides previous and current vertex lies on

			if(b_good_side == b_prev_good_side) {
				if(b_good_side) {
					// both sides in front, keep vertex in the list 'as is'
				} else {
					// both sides in back, remove current vertex from the list

					m_vertex_list.erase(m_vertex_list.begin() + i);
					-- i;
					-- n;
				}
			} else {
				float f_t = r_t_plane.f_Vector_Dist(t_vertex);
				float f_subdivide_thresh = 1.0f / (f_t - r_t_plane.f_Vector_Dist(t_prev_vertex));
				f_t *= f_subdivide_thresh;
				// calculate distance of subdivision along the edge

				f_subdivide_thresh = fabs(f_subdivide_thresh * f_epsilon_ex);
				if(f_t > f_subdivide_thresh && f_t < 1 - f_subdivide_thresh) { // this is still meaningful here
					// make sure that no tiny edges are created by the split

					TVertStruct t_split_vertex(t_vertex);
					t_split_vertex = t_vertex * (1 - f_t) + t_prev_vertex * f_t;
					// calculate a new vertex lying on edge in place where it's cut by the plane

					if(b_good_side) {
						// prev side was on the wrong position of plane, we have to insert
						// split vertex before current vertex

						m_vertex_list.insert(m_vertex_list.begin() + i, t_split_vertex);
						++ i;
						++ n;
					} else {
						// this side is on the wrong position, replace the current vertex by
						// split vertex

						m_vertex_list[i] = t_split_vertex;
					}
				} else if(!b_good_side) {
					m_vertex_list.erase(m_vertex_list.begin() + i);
					-- i;
					-- n;
				}
			}

			b_prev_good_side = b_good_side;
			t_prev_vertex = t_vertex;
		}
		// cut polygon

		// note that this leaves onplane vertices; these should not be deleted!
		// (looks awkward in 2D but is neccessary in 3D when split plane might be colinear with the polygon)
	}

	/**
	 *	@brief splits the polygon by a plane to two polyons
	 *
	 *	@param[in] r_t_plane is the splitting plane
	 *	@param[out] r_other is the other part of the split
	 *	@param[in] n_desired_half is the desired half of the polygon to remain after
	 *		cutting (plane_Front and plane_Back are the only allowed values),
	 *		the other half will go to r_other
	 *	@param[in] f_epsilon_ex is half of the plane thickness (edges shorter
	 *		than that will not be split)
	 *
	 *	@return Returns the position of the original polygon relative to the plane.
	 *
	 *	@note This function throws std::bad_alloc.
	 *	@note If only a single side of the split is required, using Cut() is faster.
	 */
	EPlanePos n_Split(const _TyPlane &r_t_plane, _TyPolygon &r_other,
		EPlanePos n_desired_half = plane_Front, float f_epsilon_ex = f_epsilon) // throw(std::bad_alloc)
	{
		_ASSERTE(n_desired_half == plane_Front || n_desired_half == plane_Back);

		r_other.m_n_material_id = m_n_material_id;
		r_other.m_n_surface_flags = m_n_surface_flags;
		// copy flags

		r_other.m_t_normal = m_t_normal;
		// copy normal

		r_other.m_vertex_list.clear();
		// erase all the vertices

		if(m_vertex_list.empty()) {
			_ASSERTE(r_other.m_vertex_list.empty()); // r_other is also empty now
			return plane_Onplane;
		}
		// nothing to do here

		EPlanePos n_other_half = (n_desired_half == plane_Front)? plane_Back : plane_Front;
		// the other ('wrong one') position

#if 1
		// this is the preferred algorithm, it is quite elegant, doesn't leave errant onplane vertices
		// the only downside is that it requires the initial sorting and that the plane classifications
		// are not reused

		TVertStruct t_prev_vertex = m_vertex_list.back();
		EPlanePos n_prev_side = r_t_plane.n_Vector_Pos(t_prev_vertex, f_epsilon_ex);

		if(n_prev_side == plane_Onplane) {
			EPlanePos n_prev_not_onplane_side;
			size_t n_prev_onplane_vertex;
			for(size_t i = m_vertex_list.size() - 1;;) {
				if(!i) {
					_ASSERTE(r_other.m_vertex_list.empty());
					return plane_Onplane;
				}
				// if all the vertices are onplane, just keep them in this and be done with it

				-- i; // begin with n - 2, we already know that n - 1 is onplane
				EPlanePos n_side;
				if((n_side = r_t_plane.n_Vector_Pos(m_vertex_list[i],
				   f_epsilon_ex)) != plane_Onplane) {
					if(!i) {
						if(n_side != n_desired_half) {
							_ASSERTE(r_other.m_vertex_list.empty());
							m_vertex_list.swap(r_other.m_vertex_list);
						}
						return n_side;
					}
					// there is only a single not-onplane vertex,
					// that one decides which side the polygon should be on

					n_prev_onplane_vertex = i + 1;
					n_prev_not_onplane_side = n_side;
					// otherwise remember the vertex

					break;
				}
				// see if there are any not-onplane vertices
			}
			// could trade this loop for storage, but then it would be a two-pass algorithm
			// this seems to be slightly faster ...

#if 1
			std::vector<TVertStruct> rot_vertex_list; // do it out of place
			rot_vertex_list.insert(rot_vertex_list.end(), m_vertex_list.begin() +
				n_prev_onplane_vertex, m_vertex_list.end());
			rot_vertex_list.insert(rot_vertex_list.end(), m_vertex_list.begin(),
				m_vertex_list.begin() + n_prev_onplane_vertex);
			rot_vertex_list.swap(m_vertex_list);
#else // 1
			size_t n = m_vertex_list.size(); // do it inplace
			m_vertex_list.insert(m_vertex_list.begin(), m_vertex_list.begin() +
				n_prev_onplane_vertex, m_vertex_list.end()); // this seems to be broken in the old STL
			// According to the C++03 ISO spec (§23.1.1, Table 67) (and §23.2.3, table 11 of
			// the C++11 ISO spec), as part of sequence requirements, the operation
			// a.insert(p, i, j) in a sequence container has this precondition:
			// i, j are not iterators into a.

			m_vertex_list.erase(m_vertex_list.begin() + n, m_vertex_list.end());
#endif // 1
			// turn the polygon so that the last vertex for sure is not onplane

			t_prev_vertex = m_vertex_list.back();
			n_prev_side = n_prev_not_onplane_side;
			_ASSERTE(n_prev_side == r_t_plane.n_Vector_Pos(t_prev_vertex, f_epsilon_ex)); // no need to measure again
			// previous vertex is now a different one!

			_ASSERTE(r_t_plane.n_Vector_Pos(m_vertex_list.front(), f_epsilon_ex) == plane_Onplane); // the first is onplane
			_ASSERTE(r_t_plane.n_Vector_Pos(m_vertex_list.back(), f_epsilon_ex) != plane_Onplane); // the last is not
			// make sure that the last one really is not onplane
		}
		// turns the polygon so that the last vertex is not onplane

		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			TVertStruct t_vertex = m_vertex_list[i];
			EPlanePos n_side = r_t_plane.n_Vector_Pos(t_vertex, f_epsilon_ex);
			// determine plane sides previous and current vertex lies on

			if(n_side != n_prev_side && n_side != plane_Onplane &&
			   n_prev_side != plane_Onplane) {
				float f_t = r_t_plane.f_Vector_Dist(t_vertex);
				float f_subdivide_thresh = 1 / (f_t - r_t_plane.f_Vector_Dist(t_prev_vertex));
				f_t *= f_subdivide_thresh;
				// calculate distance of subdivision along the edge

				// note that the threshold might not be required; edge split is always
				// longer or equal to plane distance, and if the plane distance is above
				// epsilon, the point will not be classified as onplane

				_ASSERTE(f_t > fabs(f_subdivide_thresh * f_epsilon_ex) &&
					f_t < 1 - fabs(f_subdivide_thresh * f_epsilon_ex));
				// make sure that no tiny edges are created by the split

				TVertStruct t_split_vertex(t_vertex);
				t_split_vertex = t_vertex * (1 - f_t) + t_prev_vertex * f_t;
				// calculate a new vertex lying on edge in place where it's cut by the plane

				m_vertex_list.insert(m_vertex_list.begin() + i, t_split_vertex);
				++ n;
				// just created a new onplane vertex, at the current index i

				t_vertex = t_split_vertex; // !!
				n_side = plane_Onplane;
			}
			// handle edge division

			if(n_side == plane_Onplane) {
				_ASSERTE(n_prev_side != plane_Onplane);
				// all the onplanes are skipped at once, n_prev_side can only
				// be plane_Onplane if n_side is not and vice versa

				EPlanePos n_next_side = plane_Onplane;
				size_t n_first_onplane = i;
				for(++ i; i < n; ++ i) {
					if((n_next_side = r_t_plane.n_Vector_Pos(m_vertex_list[i],
					   f_epsilon_ex)) != plane_Onplane)
						break;
					t_vertex = m_vertex_list[i];
				}
				-- i;
				_ASSERTE(r_t_plane.n_Vector_Pos(t_vertex, f_epsilon_ex) == plane_Onplane);
				_ASSERTE(n_next_side != plane_Onplane); // the last vertex is not onplane, the polygon is rotated that way at the beginning
				// skip all the onplanes, keep the last one in current vertex for the next loop iteration

				if(n_prev_side == n_desired_half) {
					// this vertex is onplane, and the previous vertex that was not onplane
					// was at the "good half": keep the onplane vertices in this polygon

					if(n_next_side == n_other_half)
						r_other.m_vertex_list.push_back(t_vertex);
					// in case the edge is actually crossing (not just touching), share it with the other polygon
				} else {
					// this vertex is onplane but the previous vertex that was not onplane
					// was not at the "good half": give the onplane vertices to the other polygon

					r_other.m_vertex_list.insert(r_other.m_vertex_list.end(),
						m_vertex_list.begin() + n_first_onplane, m_vertex_list.begin() + (i + 1));
					// send onplane vertices to the other polygon

					if(n_next_side == n_desired_half) {
						m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane,
							m_vertex_list.begin() + i); // keep the last one
						n -= i - n_first_onplane;
						_ASSERTE(n == m_vertex_list.size());
						i = n_first_onplane;
					} else {
						m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane,
							m_vertex_list.begin() + (i + 1));
						n -= (i + 1) - n_first_onplane;
						_ASSERTE(n == m_vertex_list.size());
						i = n_first_onplane - 1;
					}
					// remove the current vertex from the list
				}
			} else if(n_side != n_desired_half) {
				// not the desired half

				r_other.m_vertex_list.push_back(t_vertex);
				// send this vertex to the other polygon

				m_vertex_list.erase(m_vertex_list.begin() + i);
				-- i;
				-- n;
				// remove current vertex from the list
			} /*else {
				// the desired half is no-op
			}*/

			n_prev_side = n_side;
			t_prev_vertex = t_vertex;
			// remember the previous side and the previous vertex,
			// track from which side of the edge we come
		}
		// cuts the polygon

		return (m_vertex_list.empty())? n_other_half : // all in the other polygon?
			(r_other.m_vertex_list.empty())? n_desired_half : // all in this polygon?
			plane_Split; // a true split
#elif 0 // 1
		// this code is, in general, working but does not handle onplane points quite correctly.
		// if there are multiple onplane points, only a single polygon should have them (otherwise
		// overlaps, intersecting edges of the split polygons or other artifacts occur).

		TVertStruct t_prev_vertex = m_vertex_list.back();
		bool b_prev_good_side = r_t_plane.n_Vector_Pos(t_prev_vertex) != n_other_half;
		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			TVertStruct t_vertex = m_vertex_list[i];
			int n_side = r_t_plane.n_Vector_Pos(t_vertex);
			bool b_good_side = n_side != n_other_half;
			// determine plane sides previous and current vertex lies on

			if(b_good_side == b_prev_good_side) {
				if(b_good_side) {
					// both sides in front, keep vertex in the list 'as is'

					if(n_side == plane_Onplane)
						r_other.m_vertex_list.push_back(t_vertex);
					// add onplane vertices to the other one as well
				} else {
					// both sides in back, remove current vertex from the list

					r_other.m_vertex_list.push_back(t_vertex);
					m_vertex_list.erase(m_vertex_list.begin() + i);
					-- i;
					-- n;
				}
			} else {
				float f_t = r_t_plane.f_Vector_Dist(t_vertex);
				float f_subdivide_thresh = 1.0f / (f_t - r_t_plane.f_Vector_Dist(t_prev_vertex));
				f_t *= f_subdivide_thresh;
				// calculate distance of subdivision along the edge

				if(f_subdivide_thresh < 0)
					f_subdivide_thresh = -f_subdivide_thresh;
				f_subdivide_thresh *= f_epsilon;
				// calculate treshold distance

				if(f_t > f_subdivide_thresh &&
				   f_t < 1 - f_subdivide_thresh) {
					// new vertex is far enough from existing vertices ...
					// (prevents from creating tiny edges)

					// note these vertices won't be onplane because edge distance
					// is always greater or equal to plane distance

					TVertStruct t_split_vertex(t_vertex);
					t_split_vertex = t_vertex * (1 - f_t) + t_prev_vertex * f_t;
					// calculate a new vertex lying on edge in place where it's cut by the plane

					if(b_good_side) {
						// prev side was on the wrong position of plane, we have to insert
						// split vertex before current vertex

						m_vertex_list.insert(m_vertex_list.begin() + i, t_split_vertex);
						r_other.m_vertex_list.push_back(t_split_vertex);
						++ i;
						++ n;
					} else {
						// this side is on the wrong position, replace the current vertex by
						// split vertex

						r_other.m_vertex_list.push_back(t_split_vertex);
						r_other.m_vertex_list.push_back(t_vertex);
						m_vertex_list[i] = t_split_vertex;
					}
				} else {
					// vertices with close intersection (not necessarily onplane)

					if(b_good_side) {
						// prev side was on the wrong position of plane, we have to insert
						// current vertex (== split) to the other poly and keep it in this one

						r_other.m_vertex_list.push_back(t_vertex);
					} else {
						// prev side was onplane (if it wasn't we'd be in split section),
						// this side is back

						r_other.m_vertex_list.push_back(t_vertex);
						m_vertex_list.erase(m_vertex_list.begin() + i);
						-- i;
						-- n;
					}
				}
			}

			b_prev_good_side = b_good_side;
			t_prev_vertex = t_vertex;
		}
		// cut polygon
#elif 0 // 1
		TVertStruct t_prev_vertex = m_vertex_list.back();
		EPlanePos n_prev_side = r_t_plane.n_Vector_Pos(t_prev_vertex, f_epsilon_ex);

		EPlanePos n_prev_not_onplane_side;
		size_t n_prev_not_onplane_vertex;
		if(n_prev_side != plane_Onplane) {
			n_prev_not_onplane_vertex = m_vertex_list.size() - 1;
			n_prev_not_onplane_side = n_prev_side;
			// the last vertex is not onplane
		} else {
			for(size_t i = m_vertex_list.size() - 1;;) {
				if(!i) {
					_ASSERTE(r_other.m_vertex_list.empty());
					return;
					// if all the vertices are onplane, just keep them in this and be done with it (r_other is empty)
				}
				// i > 0

				-- i; // here
				// begin with n - 2, we already know that n - 1 is onplane

				EPlanePos n_side;
				if((n_side = r_t_plane.n_Vector_Pos(m_vertex_list[i],
				   f_epsilon_ex)) != plane_Onplane) {
					n_prev_not_onplane_vertex = i;
					n_prev_not_onplane_side = n_side;
					break;
				}
				// see if there are any not-onplane vertices
			}
			// could trade this loop for storage, but then it would be a two-pass algorithm
			// this seems to be slightly faster ...
		}
		// finds the previous vertex that is not onplane
		// (note that the algorithm would be simplified if we rotated the whole polygon so
		// that the last vertex is not an onplane one)

		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			TVertStruct t_vertex = m_vertex_list[i];
			EPlanePos n_side = r_t_plane.n_Vector_Pos(t_vertex, f_epsilon_ex);
			// determine plane sides previous and current vertex lies on

			if(n_side == n_prev_side && n_side != plane_Onplane) {
				if(n_side == n_desired_half) {
					// both sides in front, keep vertex in the list 'as is'

					// note there are no onplane vertices now
				} else {
					_ASSERTE(n_side == n_other_half); // not onplane here
					// both sides in back

					r_other.m_vertex_list.push_back(t_vertex);
					// send this vertex to the other polygon

					m_vertex_list.erase(m_vertex_list.begin() + i);
					_ASSERTE(!i || n_prev_not_onplane_vertex != i); // sometimes fails with i == 0
					if(n_prev_not_onplane_vertex >= i)
						-- n_prev_not_onplane_vertex; // sets to size_t(-1)
					-- i;
					-- n;
					// remove current vertex from the list
				}
			} else {
				if(n_side != plane_Onplane && n_prev_side != plane_Onplane) {
					float f_t = r_t_plane.f_Vector_Dist(t_vertex);
					float f_subdivide_thresh = 1 / (f_t - r_t_plane.f_Vector_Dist(t_prev_vertex));
					f_t *= f_subdivide_thresh;
					// calculate distance of subdivision along the edge

					// note that the threshold might not be required; edge split is always longer or equal to plane distance,
					// and if the plane distance is above epsilon, the point will not be classified as onplane

					_ASSERTE(f_t > fabs(f_subdivide_thresh * f_epsilon_ex) &&
						f_t < 1 - fabs(f_subdivide_thresh * f_epsilon_ex));
					// new vertex is far enough from existing vertices ...
					// (prevents from creating tiny edges)

					TVertStruct t_split_vertex(t_vertex);
					t_split_vertex = t_vertex * (1 - f_t) + t_prev_vertex * f_t;
					// calculate a new vertex lying on edge in place where it's cut by the plane

					m_vertex_list.insert(m_vertex_list.begin() + i, t_split_vertex);
					//if(n_prev_not_onplane_vertex != size_t(-1) && n_prev_not_onplane_vertex >= i)
					//	++ n_prev_not_onplane_vertex; // overwritten below
					++ n;
					// just created a new onplane vertex, at the current index i

					n_prev_not_onplane_side = n_prev_side;
					n_prev_not_onplane_vertex = (i + n - 1) % n; // previous one

					t_vertex = t_split_vertex; // !!
					n_side = plane_Onplane;
				}

				if(n_side == plane_Onplane) {
					if(n_prev_not_onplane_side == n_desired_half) {
						// this vertex is onplane, and the previous vertex that was not onplane was at the "good half"
						// keep the onplane vertices in this polygon
					} else {
						// this vertex is onplane but the previous vertex that was not onplane was not at the "good half"
						// give the onplane vertices to the other polygon

						r_other.m_vertex_list.push_back(t_vertex);
						// send this vertex to the other polygon

						m_vertex_list.erase(m_vertex_list.begin() + i);
						_ASSERTE(n_prev_not_onplane_vertex != i); // can't fail, i is onplane
						if(n_prev_not_onplane_vertex > i)
							-- n_prev_not_onplane_vertex;
						-- i;
						-- n;
						// remove current vertex from the list
					}
				} else if(n_side == n_desired_half) {
					if(n_prev_side == plane_Onplane) {
						if(n_prev_not_onplane_side == n_desired_half)
							r_other.m_vertex_list.push_back(t_prev_vertex);
						else {
							_ASSERTE(n_prev_not_onplane_vertex != size_t(-1)); // want to actually use it
							if(!r_other.m_vertex_list.empty()) {
								// the onplane vertex was already removed from this polygon and put to the other
								m_vertex_list.insert(m_vertex_list.begin() + i, r_other.m_vertex_list.back()); // "undelete" it - not very cool
								//printf("using other back\n"); // seems to always work, even if the order of vertices is reversed
							} else {
								size_t j = (n_prev_not_onplane_vertex + 1) % n;
								while(r_t_plane.n_Vector_Pos(m_vertex_list[(j + 1) % n], f_epsilon_ex) == plane_Onplane) // guaranteed to stop looping, n_side = n_desired_half != plane_Onplane
									j = (j + 1) % n;
								m_vertex_list.insert(m_vertex_list.begin() + i, m_vertex_list[j]); // this fails sometimes, as the next vertex might not be the last one onplane
								//printf("using prev not onplane\n"); // works with the fixup above
								// the onplane vertex was not processed yet, just copy it since
								// it is going to be deleted later
							}
							if(n_prev_not_onplane_vertex >= i)
								++ n_prev_not_onplane_vertex;
							++ i;
							++ n;
						}
						// send this vertex to the other polygon
					}
					// the previous vertex was onplane, this one is not: share one vertex with the other side

					// keep it in this polygon
				} else { // not the desired half
					if(n_prev_not_onplane_side == n_desired_half) {
						// previous not-onplane was on the desired half, the onplane vertices are in this

						r_other.m_vertex_list.push_back(t_prev_vertex);
						// send this vertex to the other polygon
					} else {
						// previous onplane was on the wrong half, the onplane vertices are in the other polygon
						// no need to share onplane vertices, this polygon will be empty
					}
					r_other.m_vertex_list.push_back(t_vertex);
					// send this vertex to the other polygon

					m_vertex_list.erase(m_vertex_list.begin() + i);
					_ASSERTE(!i || n_prev_not_onplane_vertex != i); // sometimes fails with i == 0
					if(n_prev_not_onplane_vertex >= i)
						-- n_prev_not_onplane_vertex; // sets to size_t(-1)
					-- i;
					-- n;
					// remove current vertex from the list
				}
			}

			if(n_side != plane_Onplane && n_prev_side == plane_Onplane) { // note that the second part of the condition should not be needed
				n_prev_not_onplane_side = n_side;
				n_prev_not_onplane_vertex = i;
			}
			n_prev_side = n_side; // !!
			t_prev_vertex = t_vertex;
		}
		// cuts the polygon

		// this implementation is not bad, except that it leaves one vertex in r_t_other even if r_t_other is empty (easy to fix, just make a counter of onplane verts in r_t_other and clear it if it equals all vertices) <- can't be fixed, don't know if the plane is actually crossed or not (can be avoided in situations when split vertex is calculated, then we know, otherwise would need n vertices lookahead)
		// also it produces different results for front-facing and back-facing polygons if there are multiple onplane points (the shared vertex is always the last one in order of vertices, not the one closer to the "good side")
		// the algorithm deletes onplane vertices too quickly, requiring putting one onplane vertex back, which is a rather complicated code <- also would work with lookahead
#elif 0 // 1
		// this produces malformed polygons

		// for each vertex: see if it is on the "good" side or onplane. record beginning of onplane
		// runs and ends of onplane runs. if one side is entirely onplane, it is deleted and the other
		// side is completely replaced by it. if there are runs of onplane vertices, only a single one
		// is shared by the polygons (splits from different sides will then yield different results!)

		TVertStruct t_prev_vertex = m_vertex_list.back();
		EPlanePos n_prev_side = r_t_plane.n_Vector_Pos(t_prev_vertex, f_epsilon_ex);

		EPlanePos n_next_not_onplane_side = plane_Onplane;
		TVertStruct t_next_nop_vertex = t_prev_vertex; // no default ctor?
		EPlanePos n_prev_not_onplane_side;
		size_t n_prev_not_onplane_vertex;
		if(n_prev_side != plane_Onplane) {
			n_prev_not_onplane_vertex = m_vertex_list.size() - 1;
			n_prev_not_onplane_side = n_prev_side;
			// the last vertex is not onplane
		} else {
			for(size_t i = m_vertex_list.size() - 1;;) {
				if(!i) {
					_ASSERTE(r_other.m_vertex_list.empty());
					return;
					// if all the vertices are onplane, just keep them in this and be done with it (r_other is empty)
				}
				// i > 0

				-- i; // here
				// begin with n - 2, we already know that n - 1 is onplane

				EPlanePos n_side;
				if((n_side = r_t_plane.n_Vector_Pos(m_vertex_list[i],
				   f_epsilon_ex)) != plane_Onplane) {
					n_prev_not_onplane_vertex = i;
					n_prev_not_onplane_side = n_side;
					break;
				}
				// see if there are any not-onplane vertices
			}
			// could trade this loop for storage, but then it would be a two-pass algorithm
			// this seems to be slightly faster ...

			for(size_t i = 1, n = m_vertex_list.size(); i < n; ++ i) {
				EPlanePos n_side;
				if((n_side = r_t_plane.n_Vector_Pos(m_vertex_list[i],
				   f_epsilon_ex)) != plane_Onplane) {
					n_next_not_onplane_side = n_side;
					t_next_nop_vertex = m_vertex_list[i - 1]; // use the last onplane
					break;
				}
				// see if there are any not-onplane vertices
			}
			// also look ahead
		}
		// finds the previous vertex that is not onplane
		// (note that the algorithm would be simplified if we rotated the whole polygon so
		// that the last vertex is not an onplane one)

		for(size_t i = 0, n = m_vertex_list.size(); i < n; ++ i) {
			TVertStruct t_vertex = m_vertex_list[i];
			EPlanePos n_side = r_t_plane.n_Vector_Pos(t_vertex, f_epsilon_ex);
			// determine plane sides previous and current vertex lies on

			if(n_side == n_prev_side && n_side != plane_Onplane) { // n_prev_side should not be plane_Onplane at all, except for the first vertex
				if(n_side == n_desired_half) {
					// both sides in front, keep vertex in the list 'as is'

					// note there are no onplane vertices now
				} else {
					_ASSERTE(n_side == n_other_half); // not onplane here
					// both sides in back

					r_other.m_vertex_list.push_back(t_vertex);
					// send this vertex to the other polygon

					m_vertex_list.erase(m_vertex_list.begin() + i);
					_ASSERTE(!i || n_prev_not_onplane_vertex != i); // sometimes fails with i == 0
					if(n_prev_not_onplane_vertex >= i)
						-- n_prev_not_onplane_vertex; // sets to size_t(-1)
					-- i;
					-- n;
					// remove the current vertex from the list
				}
			} else {
				if(n_side != plane_Onplane && n_prev_side != plane_Onplane) {
					float f_t = r_t_plane.f_Vector_Dist(t_vertex);
					float f_subdivide_thresh = 1 / (f_t - r_t_plane.f_Vector_Dist(t_prev_vertex));
					f_t *= f_subdivide_thresh;
					// calculate distance of subdivision along the edge

					// note that the threshold might not be required; edge split is always longer or equal to plane distance,
					// and if the plane distance is above epsilon, the point will not be classified as onplane

					_ASSERTE(f_t > fabs(f_subdivide_thresh * f_epsilon_ex) &&
						f_t < 1 - fabs(f_subdivide_thresh * f_epsilon_ex));
					// new vertex is far enough from existing vertices ...
					// (prevents from creating tiny edges)

					TVertStruct t_split_vertex(t_vertex);
					t_split_vertex = t_vertex * (1 - f_t) + t_prev_vertex * f_t;
					// calculate a new vertex lying on edge in place where it's cut by the plane

					m_vertex_list.insert(m_vertex_list.begin() + i, t_split_vertex);
					if(n_prev_not_onplane_vertex != size_t(-1) && n_prev_not_onplane_vertex >= i)
						++ n_prev_not_onplane_vertex;
					++ n;
					// just created a new onplane vertex, at the current index i

					t_vertex = t_split_vertex; // !!
					n_side = plane_Onplane;
				}

				if(n_side == plane_Onplane) {
					size_t n_first_onplane = i;
					EPlanePos n_from_side = (n_prev_side != plane_Onplane)? n_prev_side : n_prev_not_onplane_side;
					EPlanePos n_to_side = plane_Onplane;
					for(++ i; i < n; ++ i) {
						n_to_side = r_t_plane.n_Vector_Pos(m_vertex_list[i], f_epsilon_ex);
						if(n_to_side != plane_Onplane)
							break;
					}
					// find the next vertex that is not onplane

					_ASSERTE(n_from_side != plane_Onplane);
					if(n_from_side == n_desired_half) {
						// vertices between n_first_onplane to i are onplane
						// (excluding i itself, which is either not onplane or equal to n)
						// keep the onplane vertices in this polygon

						if(n_to_side == n_other_half) {
							r_other.m_vertex_list.push_back(m_vertex_list[i - 1]);
							// share the last onplane vertex with the other polygon
						} else if(n_to_side == plane_Onplane) {
							// the range wraps arround, need to know which side the polygon crosses to

							_ASSERTE(n_next_not_onplane_side != plane_Onplane); // this should be set in this case
							if(n_next_not_onplane_side == n_other_half)
								r_other.m_vertex_list.push_back(t_next_nop_vertex); // can't be solved differently, the vertex does not exist anymore (not sure which polygon it is in, index is unknown)
							// share the last onplane vertex with the other polygon
						} else {
							// it is coming back to n_desired_half, the plane is just being
							// touched. keep the vertex, no need to share it.
						}
					} else {
						_ASSERTE(n_from_side == n_other_half);
						// vertices between n_first_onplane to i are onplane
						// (excluding i itself, which is either not onplane or equal to n)
						// give the onplane vertices to the other polygon

						r_other.m_vertex_list.insert(r_other.m_vertex_list.end(),
							m_vertex_list.begin() + n_first_onplane, m_vertex_list.begin() + i);
						// copy the onplane vertices to the other polygon

						if(n_to_side == n_desired_half) {
							if(n_prev_side != plane_Onplane || (!n_first_onplane &&
							   n_prev_not_onplane_vertex == n - 1)) {
								m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane + 1,
									m_vertex_list.begin() + i);
								// the n_first_onplane is the first one, just do not delete it
							} else {
								// the first onplane is at the end of the polygon

								_ASSERTE(n_next_not_onplane_side == n_desired_half); // this should be set in this case

								// this vertex will be handled later
							}
							// share the last onplane vertex with this polygon
						} else if(n_to_side == plane_Onplane) {
							// range wraps arround

							_ASSERTE(n_next_not_onplane_side != plane_Onplane); // this should be set in this case

							if(n_next_not_onplane_side == n_desired_half) {
								m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane + 1,
									m_vertex_list.begin() + i);
								// the n_first_onplane is the first one, just do not delete it
							} else {
								m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane, m_vertex_list.begin() + i);
								// it is touching the plane from the other side, this polygon gets nothing
							}
						} else {
							m_vertex_list.erase(m_vertex_list.begin() + n_first_onplane, m_vertex_list.begin() + i);
							// only in the other polygon
						}
					}

					-- i; // to compensate the ++ i in the loop
					n = m_vertex_list.size();
				} else if(n_side == n_desired_half) {
					// keep it in this polygon
				} else { // not the desired half
					r_other.m_vertex_list.push_back(t_vertex);
					// send this vertex to the other polygon

					m_vertex_list.erase(m_vertex_list.begin() + i);
					-- i;
					-- n;
					// remove current vertex from the list
				}
			}

			n_prev_side = n_side; // !!
			t_prev_vertex = t_vertex;
		}
		// cuts the polygon

		// this algo requires rather complex look-back and look-ahead handling
		// either turn vertices in this polygon so that any onplane run is not crossing over the end
		// or do a split to two new polygons while keeping the old vertex list and then swap (or even output two!)

		// also, it is *not* a good idea to share the first onplane vertex with the other polygon,
		// it must *always* be *the last* one!
#endif // 1
	}

protected:
	/**
	 *	@brief performs a polygon-ray hit test
	 *
	 *	@tparam b_cull_front_facing is front-facing ray cull flag
	 *	@tparam b_cull_back_facing is back-facing ray cull flag
	 *	@tparam CIntersectInfo is type of intersection info to calculate and return
	 *
	 *	@param[in] r_v_org is ray origin
	 *	@param[in] r_v_dir is ray direction
	 *	@param[in] f_epsilon_ex is epsilon
	 *
	 *	@return Returns true on hit, false on miss.
	 */
	template <const bool b_cull_front_facing, const bool b_cull_back_facing, class CIntersectInfo>
	bool b_RayHit(const Vector3f &r_v_org, const Vector3f &r_v_dir,
		CIntersectInfo &r_t_iinfo, float f_epsilon_ex = 1e-5f) const
	{
		if(m_vertex_list.size() < 3)
			return false;
		Vector3f v_vert0 = (Vector3f)m_vertex_list.front();
		Vector3f v_vert1 = (Vector3f)m_vertex_list[1];
		for(size_t i = 2, n = m_vertex_list.size(); i < n; ++ i) {
			Vector3f v_vert2 = (Vector3f)m_vertex_list[i];

			Vector3f v_edge1 = v_vert1 - v_vert0,
					 v_edge2 = v_vert2 - v_vert0;
			v_vert1 = v_vert2; // not used below
			Vector3f v_p = r_v_dir.v_Cross(v_edge2);
			float f_det = v_edge1.f_Dot(v_p);

			if(b_cull_front_facing) {
				if(f_det < f_epsilon_ex)
					continue;
			}
			if(b_cull_back_facing) {
				if(f_det > -f_epsilon_ex)
					continue;
			}
			if(!b_cull_back_facing && !b_cull_back_facing) {
				if(f_det > -f_epsilon_ex && f_det < f_epsilon_ex)
					continue;
			}
			// branches are const, should get optimized away

			float f_inv_det = 1 / f_det;

			Vector3f v_t = r_v_org - v_vert0;

			float f_u = v_t.f_Dot(v_p) * f_inv_det;
			if(f_u < 0 || f_u > 1)
				continue;

			Vector3f v_q = v_t.v_Cross(v_edge1);

			float f_v = r_v_dir.f_Dot(v_q) * f_inv_det;
			if(f_v < 0 || f_u + f_v > 1)
				continue;

			r_t_iinfo = CIntersectInfo(v_edge2, v_q, f_inv_det, f_u, f_v, i);
			// calculate / copy any and all intersection info that is needed

			return true;
		}
		return false;
	}
};

#endif //__POLYGON_INCLUDED
