/*
 * $Id: geometrytriangles.h 4681 2017-09-05 14:35:37Z antonin.misek-skoda-auto.cz $
 *
 * Description  : One line description of file.
 * Author       : Vaclav Kyba <mail/Jabber: vaseo1@gmail.com>  <ICQ: 98576293>
 *
 * Purpose      :
 *     Long description of what the file is for.
 */


#ifndef __GEOMETRY_TRIANGLES__H__
#define __GEOMETRY_TRIANGLES__H__

//#include <vector>
#include <wx/sstream.h>
//#include "common.h"
#include "flexilog.h"
#include "geometry.h"
#include <wx/txtstrm.h>
#include <algorithm>
#include "VrutApi.h"

int VRUT_CORE_API geometrytriangles_cmp(const void* a, const void* b );

namespace VRUT
{
	///Container for triangulated geometry
	class GeometryTriangles : public Geometry
	{
	public:
		static const GEOMETRY_TYPE ClassType = GEOMETRY_TRIANGLE;

		typedef unsigned int Index;
		typedef Vector2 TexCoord;
		///Primitive organization
		enum PRIMITIVE_TYPE
		{
			TRI_LIST,	///<Triangle list
			TRI_FAN,	///<Triangle fan
			TRI_STRIP,	///<Triangle strip
			POLYGON,	///<Polygon
			LINES,		///<Lines
			LINE_STRIP,	///<Line strip
			POINTS,		///<Points
			QUADS		///<Quads
		};
		typedef std::pair<Index, PRIMITIVE_TYPE> TriDesc;
		typedef std::vector<std::pair<Index, PRIMITIVE_TYPE> > TriDescList;

		// ************************************************************************
		// **************	Start data for NURBS	*******************************
		// ************************************************************************

		// This entities are for P section

		/// Entity 100
		/// Circular arc entity
		struct CircularArcEntity
		{
			float zt,x1,y1,x2,y2,x3,y3;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 102
		/// Composite Curve
		struct CompositeCurve
		{
			std::vector<int> de;	// pointers on DE

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					Entity 106 and his forms			----
		// -------------------------------------------------------------


		struct AIEForm1		// for form number 1-3,11-13
		{
			int n,ip;
			float zt;
			std::vector<float> x,y,z,i,j,k;
		};

		struct AIEForm2		// for form number 20,21,31-38,40,63
		{
			int n,ip;
			float zt;
			std::vector<float> x,y;
		};


		///Entity 106
		struct Entity106
		{
			int formType;
			AIEForm1 form1;
			AIEForm2 form2;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 106 and his forms		----
		// -------------------------------------------------------------	

		/// Entity108
		/// Unbounded Plane Entity
		struct PlaneEntity
		{
			float a,b,c,d,x,y,z,size;
			int ptr;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 110
		/// Line Entity
		struct LineEntity
		{
			float x1,y1,z1,x2,y2,z2;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 116
		/// Point
		struct PointEntity
		{
			Vector3 point;
			int ptr;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 120
		/// Surface of revolution entity
		struct SurfaceOfRevolutionEntity
		{
			float sa,ta;
			int l,c;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 124
		/// Transformation matrix
		struct TransformationMatrix
		{
			Matrix mat;

			// additional pointers
			std::vector<int> deNa, deNp;
			
			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 126
		/// Rational B-Spline Curve
		struct RationalBSplineCurve
		{
			int k;
			int m;
			int prop1, prop2, prop3, prop4;
			std::vector<float> t;		// knot sequence
			std::vector<float> w;		// weight		
			std::vector<Vector3> controlPoints;	//control points
			float v0,v1;
			Vector3 norm;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 

			wxString ToString() const
			{
				wxString res(wxString::Format(wxT("<RationalBSplineCurve> k=%i, m=%i, v0=%f, v1=%f, t=("), k, m, v0, v1));
				res+=TextWriter::WriteToWxString(&t[0], t.size(), wxT(", "))+wxT("), w=(");
				res+=TextWriter::WriteToWxString(&w[0], w.size(), wxT(", "))+wxT("), controlPoints=(");
				for (auto &controlPoint : controlPoints)
					res+=wxT("(")+TextWriter::WriteToWxString(controlPoint._v, 3, wxT(", "))+wxT(")");
				res+=wxT(")");
				return res;
			}
		};

		/// Entity 128
		/// Rational B-Spline Surface
		struct RationalBSplineSurface
		{
			int k1,k2;
			int m1,m2;
			int prop1,prop2,prop3,prop4,prop5;
			std::vector<float> s;		// first knot sequence
			std::vector<float> t;		// second knot sequence
			std::vector<float> w;		// weight		
			std::vector<Vector3> controlPoints;	//control points
			float u0,u1,v0,v1;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section

			wxString ToString() const
			{
				wxString res(wxString::Format(wxT("<RationalBSplineSurface> k1=%i, k2=%i, m1=%i, m2=%i, u0=%f, u1=%f, v0=%f, v1=%f, t=("), k1, k2, m1,
						m2, u0, u1, v0, v1));
				res+=TextWriter::WriteToWxString(&t[0], t.size(), wxT(", "))+wxT("), s=(");
				res+=TextWriter::WriteToWxString(&s[0], s.size(), wxT(", "))+wxT("), w=(");
				res+=TextWriter::WriteToWxString(&w[0], w.size(), wxT(", "))+wxT("), controlPoints=(");
				for (auto &controlPoint : controlPoints)
					res+=wxT("(")+TextWriter::WriteToWxString(controlPoint._v, 3, wxT(", "))+wxT(")");
				res+=wxT(")");
				return res;
			}

		};

		/// Entity 142
		/// Curve on e parametric surface entity
		struct CurveOnParametricSurface
		{
			int crtn, sptr, bptr, cptr, pref;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 144
		/// Trimmed (Parametric) Surface
		struct TrimmedSurface
		{
			int pts, n1, pto;
			std::vector<int> pti;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 210
		/// General Label Entity
		struct GeneralLabelEntity
		{
			int denote,n;
			std::vector<int> dearrw;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					Entity 304 and his forms			----
		// -------------------------------------------------------------

		struct LFDEForm1
		{
			int l1;
			float l2,l3;
		};

		struct LFDEForm2
		{
			std::vector<float> l;
			std::string b;
		};

		/// Entity 304
		/// Line Font Definition Entity
		struct LineFontDefinitionEntity
		{
			int formType;
			int m;

			LFDEForm1 form1;
			LFDEForm2 form2;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 304 and his forms			----
		// -------------------------------------------------------------

		/// Entity 308
		/// Subfigure definition
		struct SubfigureDefinition
		{
			int depth;
			wxString name;
			std::vector<int> de;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		/// Entity 314
		/// Color Definition
		struct ColorDefinition
		{
			Vector3 color;
			wxString cName;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					Entity 402 and his forms			----
		// -------------------------------------------------------------

		struct AIEForm7
		{
			std::vector<int> de;
		};

		struct AIEForm3
		{
			std::vector<int> dev;
			std::vector<int> de;
		};

		///Entity402
		/// Associativity Instance Entity
		struct AssociativityInstanceEntity
		{
			int formType;
			AIEForm7 form7;
			AIEForm3 form3;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 402 and his forms		----
		// -------------------------------------------------------------	
		// -------------------------------------------------------------
		// ----					Entity 404 and his forms			----
		// -------------------------------------------------------------

		struct DrEForm0
		{			
			std::vector<int> vptr;		// size of 2
			std::vector<float> xorigin,yorigin;
			std::vector<int> dptr;	
		};

		struct DrEForm1
		{
			std::vector<int> vptr;		// size of 2
			std::vector<float> xorigin,yorigin,angle;
			std::vector<int> dptr;	
		};

		/// Entity404
		/// Drawing entity
		struct DrawingEntity
		{
			int formType;

			DrEForm0 form0;
			DrEForm1 form1;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 404 and his forms		----
		// -------------------------------------------------------------	
		// -------------------------------------------------------------
		// ----					Entity 406 and his forms			----
		// -------------------------------------------------------------
		
		// form 1
		struct PEForm1
		{
			std::vector<int> l;
		};

		///Form 12
		struct PEForm12
		{
			std::vector<wxString> names;
		};

		///Form 15
		struct PEForm15
		{
			wxString name;
		};

		///Form 16
		struct PEForm16
		{
			float xs;
			float ys;
		};

		// form 17
		struct PEForm17
		{
			int flag;
			wxString unit;
		};

		///Form 18
		struct PEForm18
		{
			float ispace;
		};

		///Entity406
		/// Property entity
		struct PropertyEntity
		{
			int formType;
			int npv;		// np in documentation

			PEForm1 form1;
			PEForm12 form12;
			PEForm15 form15;
			PEForm16 form16;
			PEForm17 form17;
			PEForm18 form18;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 406 and his forms		----
		// -------------------------------------------------------------

		///Entity 408
		/// Singular Subfigure Instance entity
		struct SingularSubfigureInstance
		{
			int formType;

			int de;
			float x;
			float y;
			float z;
			float s;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					Entity 410 and his forms		----
		// -------------------------------------------------------------

		struct VEForm0
		{
			int vno;
			float scale;
			int xvminp,yvminp,zvminp,xvmaxp,yvmaxp,zvmaxp;
		};

		struct VEForm1
		{
			int vno;
			float scale;
			
			float vpnx,vpny,vpnz;
			float vrpx,vrpy,vrpz;
			float cpx,cpy,cpz;
			float vupx,vupy,vupz;
			float vpd,umin,umax,vmin,vmax;
			int dci;
			float wmin,wmax;
		};

		///Entity 410
		/// View Entity
		struct ViewEntity
		{
			int formType;

			VEForm0 form0;
			VEForm1 form1;

			// additional pointers
			std::vector<int> deNa, deNp;

			// section d
			int dePointer;			// pointer on his DE section 
		};

		// -------------------------------------------------------------
		// ----					End entity 410 and his forms		----
		// -------------------------------------------------------------

		struct NurbsGeometry
		{
			wxString rootNodeName;			// node name i get from Entity label
			int type;
			int colorNumber;
			int transformationMatrix;

			// Circular arc entity 100
			std::vector<CircularArcEntity> e100;
			/// Composite curve (102)
			std::vector<CompositeCurve> e102;
			/// Entity (106)
			std::vector<Entity106> e106;
			/// Plane entity (108)
			std::vector<PlaneEntity> e108;
			/// Line entity (110)
			std::vector<LineEntity> e110;
			/// Point (116)
			std::vector<PointEntity> e116;
			/// Surface Of Revolution Entity (120)
			std::vector<SurfaceOfRevolutionEntity> e120;
			/// Transformation matrix (124)
			std::vector<TransformationMatrix> e124;
			/// Rational B-Spline curve (126) 
			std::vector<RationalBSplineCurve> e126;
			/// Rational B-Spline surface (128)
			std::vector<RationalBSplineSurface> e128;
			/// Curve on parametric surface (142)
			std::vector<CurveOnParametricSurface> e142;
			/// Trimmed surface (144)
			std::vector<TrimmedSurface> e144;
			/// Trimmed surface (212)
			std::vector<GeneralLabelEntity> e212;
			/// LineFontDefinitionEntity (304)
			std::vector<LineFontDefinitionEntity> e304;
			/// Subfigure definition (308)
			std::vector<SubfigureDefinition> e308;
			/// Color definition (314)
			std::vector<ColorDefinition> e314;
			/// Associativity Instance Entity (402)
			std::vector<AssociativityInstanceEntity> e402;
			/// Drawing Entity (404)
			std::vector<DrawingEntity> e404;
			/// Singular Subfigure Instance entity (408)
			std::vector<SingularSubfigureInstance> e408;
			/// property entity (406)
			std::vector<PropertyEntity> e406;
			/// View entity (410)
			std::vector<ViewEntity> e410;
			wxString ToString() const
			{
				return wxString::Format(wxT("type=%i, rootNodeName="), type)+rootNodeName;
			}
		};

		//each record is one NURBS geometry, it can consist of more elements and is independent on others
		std::vector<NurbsGeometry> nurbsGeometry;

		// ********************************************************************
		// **************	End data for NURBS	*******************************
		// ********************************************************************



		///Vertices
		std::vector<Vector3> vertices;
		///Indices
		std::vector<Index> indices;
		///Normals
		std::vector<Vector3> normals;
		///Texture coordinates
		std::vector<TexCoord> texCoords;
		///Each element defines primitive organization and index to \c Indices where it starts
		TriDescList triDescList;

		///Class constructor
		GeometryTriangles(GEOMETRY_ID _gid) : Geometry(_gid, ClassType)
		{
		}

		///Copy constructor
		GeometryTriangles(const GeometryTriangles & g) : Geometry(g)
		{			
			indices.assign(g.indices.begin(), g.indices.end());
			normals.assign(g.normals.begin(), g.normals.end());
			texCoords.assign(g.texCoords.begin(), g.texCoords.end());
			triDescList.assign(g.triDescList.begin(), g.triDescList.end());
			vertices.assign(g.vertices.begin(), g.vertices.end());
			nurbsGeometry.assign(g.nurbsGeometry.begin(), g.nurbsGeometry.end());
		}

		virtual ~GeometryTriangles() {}

		///Build AABB - Geometry overload
		virtual void BuildAABB()
		{
			if (vertices.empty())
			{
				aabb.SetInvalid();
				return;
			}
#ifdef BUILD_AABB_FROM_ALL_VERTICES
			std::vector<Vector3>::const_iterator v = vertices.begin();
			aabb.MinBound = aabb.MaxBound = *v;
			Vector3 * minB = &aabb.MinBound;
			Vector3 * maxB = &aabb.MaxBound;
			for (++v; v != vertices.end(); ++v)
			{
				if (v->x < minB->x)
					minB->x = v->x;
				else if (v->x > maxB->x)
					maxB->x = v->x;
				if (v->y < minB->y)
					minB->y = v->y;
				else if (v->y > maxB->y)
					maxB->y = v->y;
				if (v->z < minB->z)
					minB->z = v->z;
				else if (v->z > maxB->z)
					maxB->z = v->z;
			}
			///Fix zero volume AABB, TODO: check if zero volume cause real problems
			if (fabs(minB->x - maxB->x) < EPSILON)
			{
				minB->x -= EPSILON;
				maxB->x += EPSILON;
			}
			if (fabs(minB->y - maxB->y) < EPSILON)
			{
				minB->y -= EPSILON;
				maxB->y += EPSILON;
			}
			if (fabs(minB->z - maxB->z) < EPSILON)
			{
				minB->z -= EPSILON;
				maxB->z += EPSILON;
			}
#else
			aabb.SetInvalid();
			Vector3 * minB = &aabb.MinBound;
			Vector3 * maxB = &aabb.MaxBound;
			if (!indices.empty())
				*minB = *maxB = vertices[indices[0]];
			for (size_t i=1; i<indices.size(); ++i)
			{
				Vector3 &v=vertices[indices[i]];
				if (v.x < minB->x)
					minB->x = v.x;
				else if (v.x > maxB->x)
					maxB->x = v.x;
				if (v.y < minB->y)
					minB->y = v.y;
				else if (v.y > maxB->y)
					maxB->y = v.y;
				if (v.z < minB->z)
					minB->z = v.z;
				else if (v.z > maxB->z)
					maxB->z = v.z;
			}
#endif
		}

		///TODO rename to primitive instead of face
		///Gets the number of vertices in the chosen face
		///\param[in] faceNumber Index in the TriDescList vector 
		///\return \c number of vertices in the face
		size_t GetPrimitiveVertexCount(TriDescList::size_type faceNumber) const
		{
			if (faceNumber < (triDescList.size() - 1))
			{
				return triDescList[faceNumber + 1].first - triDescList[faceNumber].first;
			}
			else if (faceNumber < triDescList.size())
			{
				return indices.size() - triDescList[faceNumber].first;
			}
			return 0;
		}

		///Get geometry as list of triangles
		virtual std::vector<Triangle> Triangulate() const
		{
			size_t maxInd = indices.size();
			std::vector<Triangle> triangleList;

			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; indiceIndex+=3)
						triangleList.push_back(Triangle(vertices[indices[indiceIndex]], vertices[indices[indiceIndex+1]], vertices[indices[indiceIndex+2]]));
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						for ( ;	indiceIndex < endInd; )
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							triangleList.push_back(Triangle(last1, last2, next));
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::TRI_STRIP:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						bool sw = false;
						for ( ;	indiceIndex < endInd; sw = !sw)
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							if (sw)
								triangleList.push_back(Triangle(last1, next, last2));
							else
								triangleList.push_back(Triangle(last1, last2, next));
							last1 = last2;
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::QUADS:
					for ( ;	indiceIndex + 3 < endInd; )
					{
						Vector3 v0(vertices[indices[indiceIndex++]]);
						Vector3 v1(vertices[indices[indiceIndex++]]);
						Vector3 v2(vertices[indices[indiceIndex++]]);
						Vector3 v3(vertices[indices[indiceIndex++]]);
						triangleList.push_back(Triangle(v0, v1, v2));
						triangleList.push_back(Triangle(v0, v2, v3));
					}
					break;
				default:
					break;
				}
			}
			return triangleList;
		}

		///Get geometry as list of triangles
		virtual std::vector<TriangleData> TriangulateData() const
		{
			size_t maxInd = indices.size();
			//count number of triangles
			size_t triCount=0;
			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					triCount+=(endInd-indiceIndex)/3;
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
				case GeometryTriangles::TRI_STRIP:
					triCount+=(endInd-indiceIndex)-2;
					break;
				case GeometryTriangles::QUADS:
					triCount+=(endInd-indiceIndex)/2;
					break;
				default:
					break;
				}
			}

			//now do the triangulation
			std::vector<TriangleData> triangleList;
			triangleList.reserve(triCount);//prealloc vector
			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; indiceIndex+=3){
						size_t maxIndex = __max(indices[indiceIndex], __max(indices[indiceIndex+1], indices[indiceIndex+2]));
						triangleList.push_back(TriangleData(vertices[indices[indiceIndex]], vertices[indices[indiceIndex+1]], vertices[indices[indiceIndex+2]]));
						if (normals.size() > maxIndex)
							triangleList.back().SetNormals(normals[indices[indiceIndex]], normals[indices[indiceIndex+1]], normals[indices[indiceIndex+2]]);
						if (texCoords.size())
							triangleList.back().SetTexCoords(texCoords[indices[indiceIndex]], texCoords[indices[indiceIndex+1]], texCoords[indices[indiceIndex+2]]);
					}
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
					{
						Vector3 last1, last2;
						Vector3 lastN1, lastN2;
						Vector2 lastTC1, lastTC2;
						if (indiceIndex + 1 < endInd)
						{
							if(texCoords.size() > indices[indiceIndex]){
								lastTC1 = texCoords[indices[indiceIndex]];
							}
							if(normals.size() > indices[indiceIndex]){
								lastN1 = normals[indices[indiceIndex]];
							}
							last1 = vertices[indices[indiceIndex++]];
							if(texCoords.size() > indices[indiceIndex]){
								lastTC2 = texCoords[indices[indiceIndex]];
							}
							if(normals.size() > indices[indiceIndex]){
								lastN2 = normals[indices[indiceIndex]];
							}
							last2 = vertices[indices[indiceIndex++]];
						}
						for ( ;	indiceIndex < endInd; )
						{
							Vector2 nextTC;
							Vector3 nextN;
							if (texCoords.size() > indices[indiceIndex])
								nextTC=texCoords[indices[indiceIndex]];
							if (normals.size() > indices[indiceIndex])
								nextN=normals[indices[indiceIndex]];
							Vector3 next(vertices[indices[indiceIndex++]]);
							triangleList.push_back(TriangleData(last1, last2, next, lastN1, lastN2, nextN, lastTC1, lastTC2, nextTC));
							last2 = next;
							lastN2 = nextN;
							lastTC2 = nextTC;
						}
					}
					break;
				case GeometryTriangles::TRI_STRIP:
					{
						Vector3 last1, last2, lastN1, lastN2;
						Vector2 lastTC1, lastTC2;
						if (indiceIndex + 1 < endInd)
						{
							if (texCoords.size() > indices[indiceIndex])
								lastTC1 = texCoords[indices[indiceIndex]];
							if (normals.size() > indices[indiceIndex])
								lastN1 = normals[indices[indiceIndex]];
							last1 = vertices[indices[indiceIndex++]];
							if (texCoords.size() > indices[indiceIndex])
								lastTC2 = texCoords[indices[indiceIndex]];
							if (normals.size() > indices[indiceIndex])
								lastN2 = normals[indices[indiceIndex]];
							last2 = vertices[indices[indiceIndex++]];
						}
						bool sw = false;
						for ( ;	indiceIndex < endInd; sw = !sw)
						{
							Vector2 nextTC;
							Vector3 nextN;
							if (texCoords.size() > indices[indiceIndex])
								nextTC=texCoords[indices[indiceIndex]];
							if (normals.size() > indices[indiceIndex])
								nextN=normals[indices[indiceIndex]];
							Vector3 next(vertices[indices[indiceIndex++]]);
							if (sw)
								triangleList.push_back(TriangleData(last1, next, last2, lastN1, nextN, lastN2, lastTC1, nextTC, lastTC2));
							else
								triangleList.push_back(TriangleData(last1, last2, next, lastN1, lastN2, nextN, lastTC1, lastTC2, nextTC));
							last1 = last2;
							last2 = next;
							lastN1 = lastN2;
							lastN2 = nextN;
							lastTC1 = lastTC2;
							lastTC2 = nextTC;
						}
					}
					break;
				case GeometryTriangles::QUADS:
					for ( ;	indiceIndex + 3 < endInd; indiceIndex+=4)
					{
						size_t maxIndex = __max(__max(indices[indiceIndex], indices[indiceIndex+1]), __max(indices[indiceIndex+2], indices[indiceIndex+3]));

						triangleList.push_back(TriangleData(vertices[indices[indiceIndex]], vertices[indices[indiceIndex+1]], vertices[indices[indiceIndex+2]]));
						if (texCoords.size() > maxIndex)
							triangleList.back().SetTexCoords(texCoords[indices[indiceIndex]], texCoords[indices[indiceIndex+1]], texCoords[indices[indiceIndex+2]]);
						if (normals.size() > maxIndex)
							triangleList.back().SetNormals(normals[indices[indiceIndex]], normals[indices[indiceIndex+1]], normals[indices[indiceIndex+2]]);

						triangleList.push_back(TriangleData(vertices[indices[indiceIndex]], vertices[indices[indiceIndex+2]], vertices[indices[indiceIndex+3]]));
						if (texCoords.size() > maxIndex)
							triangleList.back().SetTexCoords(texCoords[indices[indiceIndex]], texCoords[indices[indiceIndex+2]], texCoords[indices[indiceIndex+3]]);
						if (normals.size() > maxIndex)
							triangleList.back().SetNormals(normals[indices[indiceIndex]], normals[indices[indiceIndex+2]], normals[indices[indiceIndex+3]]);
					}
					break;
				default:
					break;
				}
			}
			return triangleList;
		}

		///Get geometry as indices for TRI_LIST
		///Attention: it ignores Points, Lines and Line strips
		///If getOnlyFirstTriangle is set, only the first triangle is returned (two for quads). This parameter was introduced because of performance reasons when only the first triangle is needed.
		void Triangulate(std::vector<Index> &_indices, bool getOnlyFirstTriangle = false) const
		{
			if (indices.empty())
				return;

			size_t maxInd = indices.size();

			//preallocate _indices
			if (getOnlyFirstTriangle && !indices.empty())
				_indices.reserve(3);
			else
			{
				size_t outIndSize=0;

				for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
				{
					size_t indiceIndex = indPrimType->first;
					GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
					size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
					endInd = __min( endInd, maxInd );

					switch (indPrimType->second)
					{
					case GeometryTriangles::TRI_LIST:
						outIndSize+=endInd-indiceIndex;
						break;
					case GeometryTriangles::TRI_FAN:
					case GeometryTriangles::POLYGON:
					case GeometryTriangles::TRI_STRIP:
						outIndSize+=(endInd-indiceIndex-2)*3;
						break;
					case GeometryTriangles::QUADS:
						outIndSize+=((endInd-indiceIndex)/4)*6;
						break;
					default:
						break;
					}
				}
				_indices.reserve(outIndSize);
			}

			//triangulate
			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; )
					{
						_indices.push_back(indices[indiceIndex++]);
						_indices.push_back(indices[indiceIndex++]);
						_indices.push_back(indices[indiceIndex++]);
					}
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
					{
						Index last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = indices[indiceIndex++];
							last2 = indices[indiceIndex++];
						}
						for ( ;	indiceIndex < endInd; )
						{
							Index next = indices[indiceIndex++];
							_indices.push_back(last1);
							_indices.push_back(last2);
							_indices.push_back(next);
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::TRI_STRIP:
					{
						Index last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = indices[indiceIndex++];
							last2 = indices[indiceIndex++];
						}
						bool sw = false;
						for ( ;	indiceIndex < endInd; sw = !sw)
						{
							Index next = indices[indiceIndex++];
							if (sw)
							{
								_indices.push_back(last1);
								_indices.push_back(next);
								_indices.push_back(last2);
							}
							else
							{
								_indices.push_back(last1);
								_indices.push_back(last2);
								_indices.push_back(next);
							}
							last1 = last2;
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::QUADS:
					for ( ;	indiceIndex + 3 < endInd; )
					{
						Index i0 = indices[indiceIndex++];
						Index i1 = indices[indiceIndex++];
						Index i2 = indices[indiceIndex++];
						Index i3 = indices[indiceIndex++];
						_indices.push_back(i0);
						_indices.push_back(i1);
						_indices.push_back(i2);
						_indices.push_back(i0);
						_indices.push_back(i2);
						_indices.push_back(i3);
					}
					break;
				default:
					break;
				}
				if (getOnlyFirstTriangle && (_indices.size() >= 3))
					break;
			}
		}

		///Flip order of points in primitives (ie. CW-CCW), it don't change points, lines and line strips and currently cannot handle some other types
		void FlipCwCcw()
		{
			size_t maxInd = indices.size();

			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; indiceIndex+=3)
					{
						Index tmp=indices[indiceIndex];
						indices[indiceIndex]=indices[indiceIndex+1];
						indices[indiceIndex+1]=tmp;
					}
					break;
				/*case GeometryTriangles::TRI_FAN://TODO: perhaps it is not possible
				case GeometryTriangles::TRI_STRIP://TODO: perhaps it is not possible
				case GeometryTriangles::POLYGON:////TODO: perhaps it is not possible
					{
						for (size_t i=0; i < (endInd-indiceIndex)/2; ++i)
						{
							Index tmp=indices[indiceIndex+i];
							indices[indiceIndex+i]=indices[endInd-1-i];
							indices[endInd-1-i]=tmp;
						}
					}
					break;*/
				case GeometryTriangles::LINES:
				case GeometryTriangles::LINE_STRIP:
				case GeometryTriangles::POINTS:
					break;
				case GeometryTriangles::QUADS://TODO: Attention: it is untested
					for ( ;	indiceIndex + 3 < endInd; indiceIndex+=4)
					{
						Index tmp=indices[indiceIndex];
						indices[indiceIndex]=indices[indiceIndex+3];
						indices[indiceIndex+3]=tmp;
						tmp=indices[indiceIndex+1];
						indices[indiceIndex+1]=indices[indiceIndex+2];
						indices[indiceIndex+2]=tmp;
					}
					break;
				default:
					LOGERROR(wxString::Format(wxT("<GeometryTriangles>FlipCwCcw not supported primitive type: %i"), indPrimType->second));
					break;
				}
			}
		}

		///Check if ray intersects geometry, get intersection point distance from ray origin if it does - Geometry overload
		///NOTE: everything must be in one space (model space by default), \c dist will be in that space too!
		virtual bool Intersects(const Ray & ray, float * dist = (float *)NULL) const
		{
			float bestDist = FLT_MAX;
			bool found = false;
			size_t maxInd = indices.size();

			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; )
					{
						Vector3 v0(vertices[indices[indiceIndex++]]);
						Vector3 v1(vertices[indices[indiceIndex++]]);
						Vector3 v2(vertices[indices[indiceIndex++]]);
						float actDist;
						if (ray.IntersectsTri(v0, v1, v2, &actDist) && actDist < bestDist)
						{
							found = true;
							bestDist = actDist;
						}
					}
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						for ( ;	indiceIndex < endInd; )
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							float actDist;
							if (ray.IntersectsTri(last1, last2, next, &actDist) && actDist < bestDist)
							{
								found = true;
								bestDist = actDist;
							}
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::TRI_STRIP:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						bool sw = false;
						for ( ;	indiceIndex < endInd; sw = !sw)
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							float actDist;
							if (sw)
							{
								if (ray.IntersectsTri(last1, next, last2, &actDist) && actDist < bestDist)
								{
									found = true;
									bestDist = actDist;
								}
							}
							else
							{
								if (ray.IntersectsTri(last1, last2, next, &actDist) && actDist < bestDist)
								{
									found = true;
									bestDist = actDist;
								}
							}
							last1 = last2;
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::QUADS:
					for ( ;	indiceIndex + 3 < endInd; )
					{
						Vector3 v0(vertices[indices[indiceIndex++]]);
						Vector3 v1(vertices[indices[indiceIndex++]]);
						Vector3 v2(vertices[indices[indiceIndex++]]);
						Vector3 v3(vertices[indices[indiceIndex++]]);
						float actDist;
						if (ray.IntersectsTri(v0, v1, v2, &actDist) && actDist < bestDist)
						{
							found = true;
							bestDist = actDist;
						}
						if (ray.IntersectsTri(v0, v2, v3, &actDist) && actDist < bestDist)
						{
							found = true;
							bestDist = actDist;
						}
					}
					break;
				default:
					break;
				}
			}

			if (found && dist)
				*dist = bestDist;

			return found;
		}

		virtual bool CastRay(const Ray & ray, RayIntersectionInfo &info) const
		{
			float bestDist = FLT_MAX;
			size_t bestIndex = (size_t)~(size_t(0));
			GeometryTriangles::TriDescList::size_type bestPrimitive;
			size_t primtitiveIndex = 0;
			size_t maxInd = indices.size();
			size_t triangleIndex=(size_t)-1;

			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType, ++primtitiveIndex)
			{
				size_t indiceIndex = indPrimType->first;
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ((endIndPrimType == triDescList.end()) ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				switch (indPrimType->second)
				{
				case GeometryTriangles::TRI_LIST:
					for ( ;	indiceIndex + 2 < endInd; )
					{
						Vector3 v0(vertices[indices[indiceIndex++]]);
						Vector3 v1(vertices[indices[indiceIndex++]]);
						Vector3 v2(vertices[indices[indiceIndex++]]);
						++triangleIndex;
						float actDist;
						if (ray.IntersectsTri(v0, v1, v2, &actDist) && actDist < bestDist)
						{
							bestIndex = triangleIndex;
							bestPrimitive = primtitiveIndex;
							info.normal = Normalize(CrossProd(v0-v1, v2-v1));
							bestDist = actDist;
						}
					}
					break;
				case GeometryTriangles::TRI_FAN:
				case GeometryTriangles::POLYGON:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						for ( ;	indiceIndex < endInd; )
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							++triangleIndex;
							float actDist;
							if (ray.IntersectsTri(last1, last2, next, &actDist) && actDist < bestDist)
							{
								bestIndex = triangleIndex;
								bestPrimitive = primtitiveIndex;
								info.normal = Normalize(CrossProd(last1-last2, next-last2));
								bestDist = actDist;
							}
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::TRI_STRIP:
					{
						Vector3 last1, last2;
						if (indiceIndex + 1 < endInd)
						{
							last1 = vertices[indices[indiceIndex++]];
							last2 = vertices[indices[indiceIndex++]];
						}
						bool sw = false;
						for ( ;	indiceIndex < endInd; sw = !sw)
						{
							Vector3 next(vertices[indices[indiceIndex++]]);
							++triangleIndex;
							float actDist;
							if (sw)
							{
								if (ray.IntersectsTri(last1, next, last2, &actDist) && actDist < bestDist)
								{
									bestIndex = triangleIndex;
									bestPrimitive = primtitiveIndex;
									info.normal = Normalize(CrossProd(last1-next, last2-next));
									bestDist = actDist;
								}
							}
							else
							{
								if (ray.IntersectsTri(last1, last2, next, &actDist) && actDist < bestDist)
								{
									bestIndex = triangleIndex;
									bestPrimitive = primtitiveIndex;
									info.normal = Normalize(CrossProd(last1-last2, next-last2));
									bestDist = actDist;
								}
							}
							last1 = last2;
							last2 = next;
						}
					}
					break;
				case GeometryTriangles::QUADS:
					for ( ;	indiceIndex + 3 < endInd; )
					{
						Vector3 v0(vertices[indices[indiceIndex++]]);
						Vector3 v1(vertices[indices[indiceIndex++]]);
						Vector3 v2(vertices[indices[indiceIndex++]]);
						Vector3 v3(vertices[indices[indiceIndex++]]);
						++triangleIndex;
						float actDist;
						if (ray.IntersectsTri(v0, v1, v2, &actDist) && actDist < bestDist)
						{
							bestIndex = triangleIndex;
							bestPrimitive = primtitiveIndex;
							info.normal = Normalize(CrossProd(v0-v1, v2-v1));
							bestDist = actDist;
						}
						++triangleIndex;
						if (ray.IntersectsTri(v0, v2, v3, &actDist) && actDist < bestDist)
						{
							bestIndex = triangleIndex;
							bestPrimitive = primtitiveIndex;
							info.normal = Normalize(CrossProd(v0-v2, v3-v2));
							bestDist = actDist;
						}
					}
					break;
				default:
					break;
				}
			}

			if (bestIndex!=(size_t)~(size_t(0))) {
				info.primitiveIndex = bestPrimitive;
				info.dist = bestDist;
				info.index = bestIndex;
				info.geometryID = GetID();
				return true;
			}
			return false;
		}

		///Add quad perpendicular to z axis
		virtual void AddUnitQuad()
		{

			GeometryTriangles::Index startInd = GeometryTriangles::Index(vertices.size());

			std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE> indType(startInd, GeometryTriangles::TRI_FAN);

			triDescList.push_back(indType);

			normals.push_back(Vector3(0,0,1));
			normals.push_back(Vector3(0,0,1));
			normals.push_back(Vector3(0,0,1));
			normals.push_back(Vector3(0,0,1));

			vertices.push_back(Vector3(-1,-1,0));
			vertices.push_back(Vector3(1,-1,0));
			vertices.push_back(Vector3(1,1,0));
			vertices.push_back(Vector3(-1,1,0));

			texCoords.push_back(GeometryTriangles::TexCoord(0,0));
			texCoords.push_back(GeometryTriangles::TexCoord(1,0));
			texCoords.push_back(GeometryTriangles::TexCoord(1,1));
			texCoords.push_back(GeometryTriangles::TexCoord(0,1));

			indices.push_back(startInd++);
			indices.push_back(startInd++);
			indices.push_back(startInd++);
			indices.push_back(startInd++);
			
			
		}

		///Add box parallel to the XYZ axes
		virtual void AddBox(const Vector3 &position, const Vector3 &size)
		{
			float box[6];
			box[0] = position._v[0];
			box[1] = position._v[1];
			box[2] = position._v[2];
			box[3] = position._v[0]+size._v[0];
			box[4] = position._v[1]+size._v[1];
			box[5] = position._v[2]+size._v[2];
			GeometryTriangles::Index startInd = GeometryTriangles::Index(vertices.size());

			std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE> indType(startInd, GeometryTriangles::QUADS);
			triDescList.push_back(indType);

			float norms[]={0.0f,0.0f,1.0f, 0.0f,0.0f,-1.0f, -1.0f,0.0f,0.0f, 1.0f,0.0f,0.0f, 0.0f,1.0f,0.0f, 0.0f,-1.0f,0.0f};
			float textC[]={0.0f,0.0f, 1.0f,0.0f, 1.0f,1.0f, 0.0f,1.0f};
			int verI[]={0,1,5, 3,1,5, 3,4,5, 0,4,5,  3,1,2, 0,1,2, 0,4,2, 3,4,2,  0,1,2, 0,1,5, 0,4,5, 0,4,2,  3,1,5, 3,1,2, 3,4,2, 3,4,5,  0,4,5, 3,4,5, 3,4,2, 0,4,2,  0,1,2, 3,1,2, 3,1,5, 0,1,5}; // 14.3.2010 FIX: probable mistake in vertex order, now it is CCW
			size_t k=0;
			for (size_t i=0; i<6; i++)
			{
				for (size_t j=0; j<4; j++)
				{
					normals.push_back(Vector3(norms+i*3));
					texCoords.push_back(GeometryTriangles::TexCoord(textC+j*2));
					vertices.push_back(Vector3(box[verI[k]], box[verI[k+1]], box[verI[k+2]]));
					k+=3;
					indices.push_back(startInd++);
				}
			}
		}

		///Add content of another geometry
		virtual bool AddGeometry(const Geometry *geom)
		{
			if (geom==NULL)
				return true;
			if (geom->type!=Geometry::GEOMETRY_TRIANGLE)
			{
				LOGERROR(wxString::Format(wxT("Cannot add geometry type %i to GEOMETRY_TRIANGLE."), geom->type));
				return false;
			}
			GeometryTriangles *geoTri = (GeometryTriangles *)geom;
			size_t vertCount=vertices.size();
			size_t indCount=indices.size();
			vertices.insert(vertices.end(), geoTri->vertices.begin(), geoTri->vertices.end());

			if (!geoTri->normals.empty())
			{
				normals.resize(vertCount, Vector3());//TODO: solve missing normals properly
				normals.insert(normals.end(), geoTri->normals.begin(), geoTri->normals.end());
			}
			if (!normals.empty())
				normals.resize(vertices.size(), Vector3());//TODO: solve missing normals properly

			if (!geoTri->texCoords.empty())
			{
				texCoords.resize(vertCount, Vector2());//TODO: solve missing texCoords properly
				texCoords.insert(texCoords.end(), geoTri->texCoords.begin(), geoTri->texCoords.end());
			}
			if (!texCoords.empty())
				texCoords.resize(vertices.size(), Vector2());//TODO: solve missing texCoords properly

			if (vertCount==0)
				indices.insert(indices.begin(), geoTri->indices.begin(), geoTri->indices.end());
			else
			{
				indices.reserve(indCount+geoTri->indices.size());
				for (auto index : geoTri->indices)
					indices.push_back((NODE_ID)(vertCount+index));
			}
			if (indCount==0)
				triDescList.insert(triDescList.begin(), geoTri->triDescList.begin(), geoTri->triDescList.end());
			else
			{
				triDescList.reserve(triDescList.size()+geoTri->triDescList.size());
				for (auto &triDesc : geoTri->triDescList)
				{
					std::pair<Index, PRIMITIVE_TYPE> indType(triDesc);
					indType.first+=(Index)indCount;
					triDescList.push_back(indType);
				}
			}
			if (!geoTri->nurbsGeometry.empty())
				nurbsGeometry.insert(nurbsGeometry.end(), geoTri->nurbsGeometry.begin(), geoTri->nurbsGeometry.end());
			return true;
		}

		virtual bool AddGeometries(const std::vector<Geometry *> &geoms)
		{
			//compute needed sizes of arrays
			size_t vertCount=vertices.size();
			bool norms=!normals.empty();
			bool texCs=!texCoords.empty();
			size_t indCount=indices.size();
			size_t triDescCount=triDescList.size();
			size_t nurbsCount=nurbsGeometry.size();

			for (auto &geom : geoms)
			{
				if (!geom)
					continue;
				if (geom->type!=Geometry::GEOMETRY_TRIANGLE)
				{
					LOGERROR(wxString::Format(wxT("Cannot add geometry type %i to GEOMETRY_TRIANGLE."), geom->type));
					continue;
				}
				GeometryTriangles *geoTri = (GeometryTriangles *)geom;
				vertCount+=geoTri->vertices.size();
				if (!norms)
					norms=!geoTri->normals.empty();
				if (!texCs)
					texCs=!geoTri->texCoords.empty();
				indCount+=geoTri->indices.size();
				triDescCount+=geoTri->triDescList.size();
				nurbsCount+=geoTri->nurbsGeometry.size();
			}

			//preallocate arrays
			vertices.reserve(vertCount);
			if (norms)
				normals.reserve(vertCount);
			if (texCs)
				texCoords.reserve(vertCount);
			indices.reserve(indCount);
			triDescList.reserve(triDescCount);
			nurbsGeometry.reserve(nurbsCount);

			//add it really
			for (auto &geom : geoms)
			{
				if (!geom || geom->type!=Geometry::GEOMETRY_TRIANGLE)
					continue;
				GeometryTriangles *geoTri = (GeometryTriangles *)geom;
				size_t actVertCount=vertices.size();
				if (norms)
				{
					normals.resize(actVertCount, Vector3());
					if (!geoTri->normals.empty())
						normals.insert(normals.end(), geoTri->normals.begin(), geoTri->normals.end());
				}
				if (texCs)
				{
					texCoords.resize(actVertCount, Vector2());
					if (!geoTri->texCoords.empty())
						texCoords.insert(texCoords.end(), geoTri->texCoords.begin(), geoTri->texCoords.end());
				}
				size_t actIndCount=indices.size();
				if (actVertCount==0)
					indices.assign(geoTri->indices.begin(), geoTri->indices.end());
				else
				{
					for (size_t j=0; j<geoTri->indices.size(); ++j)
						indices.push_back((NODE_ID)(actVertCount+geoTri->indices.at(j)));
				}
				if (actIndCount==0)
					triDescList.assign(geoTri->triDescList.begin(), geoTri->triDescList.end());
				else
				{
					for (size_t j=0; j<geoTri->triDescList.size(); ++j)
					{
						std::pair<Index, PRIMITIVE_TYPE> indType(geoTri->triDescList.at(j));
						indType.first+=(Index)actIndCount;
						triDescList.push_back(indType);
					}
				}
				vertices.insert(vertices.end(), geoTri->vertices.begin(), geoTri->vertices.end());
				if (!geoTri->nurbsGeometry.empty())
					nurbsGeometry.insert(nurbsGeometry.end(), geoTri->nurbsGeometry.begin(), geoTri->nurbsGeometry.end());
			}
			size_t actVertCount=vertices.size();
			if (norms)
				normals.resize(actVertCount, Vector3());
			if (texCs)
				texCoords.resize(actVertCount, Vector2());
			return true;
		}

		//Transform vertices and flush transformation matrix
		virtual bool Transform(const Matrix *matrix)
		{
			for (auto &vertex : vertices)
				vertex = matrix->TransformCoordHom(vertex);
			for (auto &normal : normals)
				normal = matrix->TransformNormal(normal);
			return true;
		}

		///method for computing normals
		///\param[in] creaseAngle is an angle to indicate maximal angle which should be considered as smooth
		//TODO: the method is not working properly at all, it must be rewriten
		void ComputeNormals(const float & creaseAngle)
		{
			normals.clear();
			float x = 0, y = 0, z = 0;
			Vector3 normal;
/*
			size_t it_vert=0;
			std::vector<Index> triInd;
			Triangulate(triInd);
			std::pair<Vector3,size_t> *verticestmp = new std::pair<Vector3, size_t>[triInd.size()];
			for (size_t i=0; i<triInd.size()-2; i+=3)
			{
				x = y = z = 0;
				for ( unsigned loop = 0; loop < 2; loop++){
					x += (vertices[triInd[it_vert]]._v[1] - vertices[triInd[it_vert+1]]._v[1]) * (vertices[triInd[it_vert]]._v[2] + vertices[triInd[it_vert+1]]._v[2]);
					y += (vertices[triInd[it_vert]]._v[2] - vertices[triInd[it_vert+1]]._v[2]) * (vertices[triInd[it_vert]]._v[0] + vertices[triInd[it_vert+1]]._v[0]);
					z += (vertices[triInd[it_vert]]._v[0] - vertices[triInd[it_vert+1]]._v[0]) * (vertices[triInd[it_vert]]._v[1] + vertices[triInd[it_vert+1]]._v[1]);
					it_vert++;
				}
				x += (vertices[triInd[it_vert]]._v[1] - vertices[triInd[it_vert - 2]]._v[1]) * (vertices[triInd[it_vert]]._v[2] + vertices[triInd[it_vert - 2]]._v[2]);
				y += (vertices[triInd[it_vert]]._v[2] - vertices[triInd[it_vert - 2]]._v[2]) * (vertices[triInd[it_vert]]._v[0] + vertices[triInd[it_vert - 2]]._v[0]);
				z += (vertices[triInd[it_vert]]._v[0] - vertices[triInd[it_vert - 2]]._v[0]) * (vertices[triInd[it_vert]]._v[1] + vertices[triInd[it_vert - 2]]._v[1]);
				//copy face normal to all vertices
#define PRECISION 0.0001
				normal = Vector3((abs(x)<PRECISION)?0:x,(abs(y)<PRECISION)?0:y,(abs(z)<PRECISION)?0:z);
				normal = normal.Normalize();
				it_vert -= 2;

				for ( unsigned j = 0; j != 3; j++){
					verticestmp[it_vert] = std::pair<Vector3,size_t> (vertices[triInd[it_vert]],it_vert);
					normals.push_back(normal);
					it_vert++;
				}
			}*/
			
			std::pair<Vector3,size_t> *verticestmp = new std::pair<Vector3, size_t>[indices.size()];
			for (GeometryTriangles::TriDescList::iterator it = triDescList.begin(); it != triDescList.end(); ++it)
			{
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = it + 1;
				size_t endInd = ( endIndPrimType == triDescList.end()) ? indices.size() : endIndPrimType->first;
				endInd = __min( endInd,  indices.size() );
				size_t it_vert = it->first;
				unsigned int count = 0;

				while (it_vert + count < endInd)
				{
					switch (it->second)
					{
					case GeometryTriangles::TRI_LIST:
						count = 2;
						break;
					case GeometryTriangles::TRI_FAN:
					case GeometryTriangles::TRI_STRIP:
					case GeometryTriangles::POLYGON:
						if ( it+1 == triDescList.end())
							count = 2;
						else count = (it+1)->first - it->first - 1;
						break;
					case LINES:
						count = 1;
						break;
					case LINE_STRIP:
					case POINTS:
						count = 0;
						break;
					case GeometryTriangles::QUADS:
						count = 3;
						break;
					default:
						LOGERROR(wxString::Format(wxT("Unsupported primitive type %i"), it->second));
						delete[] verticestmp;
						return;
					}

					//Newell's method for computing surface face normals
					//it is a robust method - dealing correctly also with nonconvex polygon
					x = y = z = 0;
					for ( unsigned loop = 0; loop < count; loop++){
						x += (vertices[indices[it_vert]]._v[1] - vertices[indices[it_vert+1]]._v[1]) * (vertices[indices[it_vert]]._v[2] + vertices[indices[it_vert+1]]._v[2]);
						y += (vertices[indices[it_vert]]._v[2] - vertices[indices[it_vert+1]]._v[2]) * (vertices[indices[it_vert]]._v[0] + vertices[indices[it_vert+1]]._v[0]);
						z += (vertices[indices[it_vert]]._v[0] - vertices[indices[it_vert+1]]._v[0]) * (vertices[indices[it_vert]]._v[1] + vertices[indices[it_vert+1]]._v[1]);
						it_vert++;
					}
					x += (vertices[indices[it_vert]]._v[1] - vertices[indices[it_vert - count]]._v[1]) * (vertices[indices[it_vert]]._v[2] + vertices[indices[it_vert - count]]._v[2]);
					y += (vertices[indices[it_vert]]._v[2] - vertices[indices[it_vert - count]]._v[2]) * (vertices[indices[it_vert]]._v[0] + vertices[indices[it_vert - count]]._v[0]);
					z += (vertices[indices[it_vert]]._v[0] - vertices[indices[it_vert - count]]._v[0]) * (vertices[indices[it_vert]]._v[1] + vertices[indices[it_vert - count]]._v[1]);
					//copy face normal to all vertices
#define PRECISION 0.0001f
					normal = Vector3((fabs(x)<PRECISION)?0:x,(fabs(y)<PRECISION)?0:y,(fabs(z)<PRECISION)?0:z);
					normal = normal.Normalize();
					it_vert -= count;

					for ( unsigned i = 0; i != count + 1; i++){
						verticestmp[it_vert] = std::pair<Vector3,size_t> (vertices[indices[it_vert]],it_vert);
						normals.push_back(normal);
						it_vert++;
					}
				}
			}

			//sort vertices, so that the same vertices are together
			qsort(verticestmp, indices.size()-1, sizeof(std::pair<Vector3,size_t>), geometrytriangles_cmp);
			//average normals for the same vertices and less angle than creaseAngle
			float cosAngle;
			Vector3 temp(0,0,0);
			std::vector<size_t> index;
			for ( size_t i = 0; i < indices.size(); ){
				temp = normals.at(verticestmp[i].second);
				size_t number = i;
				while ((number+1<indices.size()) && (verticestmp[number].first == verticestmp[number+1].first)){
					cosAngle = normals.at(verticestmp[number].second).Dot(normals.at(verticestmp[number+1].second));
					cosAngle = cosAngle / (normals.at(verticestmp[number].second).Length() * normals.at(verticestmp[number+1].second).Length());
					if ( creaseAngle > acos(cosAngle)){
						index.push_back(verticestmp[number].second);
						index.push_back(verticestmp[number+1].second);
						temp += normals.at(verticestmp[number+1].second);
					}
					++number;
				}
				//copy the average normal
				if ( i < number){
					temp = temp.Normalize();
					for (auto ind : index)
						normals[ind] = temp;
					index.clear();
					i = number;
				}
				else
					++i;
			}
			delete [] verticestmp;
		}

		///Equal operator, TODO: implementation needs to be updated to handle the same geometries with different order of vertices or indices or primitives or different primitives with the same result...
		virtual bool IsSameAs(const smart_pointer<Geometry> & geo) const
		{
			if (geo->type != Geometry::GEOMETRY_TRIANGLE)
				return false;
			GeometryTriangles * geoTri = (GeometryTriangles*)geo.getPtr();
			if ((vertices.size()!=geoTri->vertices.size()) || (normals.size()!=geoTri->normals.size())
				|| (texCoords.size()!=geoTri->texCoords.size()) || (triDescList.size()!=geoTri->triDescList.size()))
				return false;
			for (size_t i=0; i<vertices.size(); i++)
				if (vertices.at(i)!=geoTri->vertices.at(i))
					return false;
			for (size_t i=0; i<normals.size(); i++)
				if (normals.at(i)!=geoTri->normals.at(i))
					return false;
			for (size_t i=0; i<texCoords.size(); i++)
				if (texCoords.at(i)!=geoTri->texCoords.at(i))
					return false;
			for (size_t i=0; i<triDescList.size(); i++)
				if ((triDescList.at(i).first!=geoTri->triDescList.at(i).first) || (triDescList.at(i).second!=geoTri->triDescList.at(i).second))
					return false;
			return true;
		}

		///Try to divide geometry by discontinuity, if whole geometry is one segment then it returns empty list
		virtual void Divide(std::vector<Geometry*> &segments) const
		{
			segments.clear();
			if (vertices.empty())
				return;
			std::vector<size_t> verticeSegment(vertices.size(), (size_t)-1);
			std::vector<std::vector<size_t> > segmentVertices(indices.size());
			size_t currentSegmentID=0, maxSegmentID=0;
			for (size_t i=0; i<triDescList.size(); ++i)
			{
				size_t startIndex=triDescList[i].first;
				size_t endIndex = (i+1<triDescList.size())?triDescList[i+1].first:indices.size();
				size_t step=endIndex-startIndex;
				switch (triDescList[i].second)
				{
				case GeometryTriangles::TRI_LIST:
					step=3;
					break;
				case GeometryTriangles::LINES:
					step=2;
					break;
				case GeometryTriangles::POINTS:
					step=1;
					break;
				case GeometryTriangles::QUADS:
					step=4;
					break;
				default:
					break;
				}

				for (size_t j=startIndex; j<endIndex; ++j)
				{
					if (j-startIndex>=step)
					{
						++maxSegmentID;
						currentSegmentID=maxSegmentID;
						startIndex=j;
					}
					if ((verticeSegment[indices[j]]!=(size_t)-1) && (verticeSegment[indices[j]]!=currentSegmentID))
					{
						//merge two segments together
						size_t newSegmentID=__min(verticeSegment[indices[j]], currentSegmentID);
						currentSegmentID=__max(verticeSegment[indices[j]], currentSegmentID);
						if (segmentVertices[newSegmentID].size()<segmentVertices[currentSegmentID].size())
						{
							size_t tmp=newSegmentID;
							newSegmentID=currentSegmentID;
							currentSegmentID=tmp;
						}
						//reassign vertices
						for (auto &segmentVertice : segmentVertices[currentSegmentID])
							verticeSegment[segmentVertice]=newSegmentID;
						segmentVertices[newSegmentID].insert(segmentVertices[newSegmentID].end(), segmentVertices[currentSegmentID].begin(), segmentVertices[currentSegmentID].end());
						segmentVertices[currentSegmentID].clear();
						//continue with newSegmentID
						currentSegmentID=newSegmentID;
						//--maxSegmentID;
					}
					else
					{
						verticeSegment[indices[j]]=currentSegmentID;
						segmentVertices[currentSegmentID].push_back(indices[j]);
					}
				}
				++maxSegmentID;
				currentSegmentID=maxSegmentID;
			}
			//check if more segments exists
			bool moreSegments=false;
			for (auto vertice : verticeSegment)
			{
				if (vertice!=0)
				{
					moreSegments=true;
					break;
				}
			}
			if (!moreSegments)
				return;
			//change segment ID to sequence 0, 1, 2, ...
			std::vector<size_t> segmentMapping(maxSegmentID, (size_t)-1);
			size_t segmentID=0;
			for (auto &vertice : verticeSegment)
			{
				if (segmentMapping[vertice]==(size_t)-1)
				{
					segmentMapping[vertice]=segmentID;
					++segmentID;
				}
				vertice=segmentMapping[vertice];
			}
			//prepare empty segments
			for (size_t i=0; i<segmentID; ++i)
				segments.push_back(new GeometryTriangles(GEOMETRY_ID_NONE));
			//create geometry for segments
			std::vector<size_t> verticeMapping(vertices.size(), (size_t)-1);
			for (size_t i=0; i<triDescList.size(); ++i)
			{
				size_t startIndex=triDescList[i].first;
				size_t endIndex = (i+1<triDescList.size())?triDescList[i+1].first:indices.size();
				size_t step=endIndex-startIndex;
				bool checkPrevType=false;
				switch (triDescList[i].second)
				{
				case GeometryTriangles::TRI_LIST:
					step=3;
					checkPrevType=true;
					break;
				case GeometryTriangles::LINES:
					step=2;
					checkPrevType=true;
					break;
				case GeometryTriangles::POINTS:
					step=1;
					checkPrevType=true;
					break;
				case GeometryTriangles::QUADS:
					step=4;
					checkPrevType=true;
					break;
				default:
					break;
				}

				for (size_t j=startIndex; j<endIndex; ++j)
				{
					if (j-startIndex>=step)
						startIndex=j;
					if (j==startIndex)
					{
						currentSegmentID=verticeSegment[indices[j]];
						if (!checkPrevType || (((GeometryTriangles*)segments[currentSegmentID])->triDescList.empty()) || (((GeometryTriangles*)segments[currentSegmentID])->triDescList.back().second!=triDescList[i].second))
							((GeometryTriangles*)segments[currentSegmentID])->triDescList.push_back(std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE> ((GeometryTriangles::Index)((GeometryTriangles*)segments[currentSegmentID])->indices.size(), triDescList[i].second));
					}
					verticeMapping[indices[j]]=((GeometryTriangles*)segments[currentSegmentID])->vertices.size();
					((GeometryTriangles*)segments[currentSegmentID])->vertices.push_back(vertices[indices[j]]);
					if (normals.size()>indices[j])
						((GeometryTriangles*)segments[currentSegmentID])->normals.push_back(normals[indices[j]]);
					if (texCoords.size()>indices[j])
						((GeometryTriangles*)segments[currentSegmentID])->texCoords.push_back(texCoords[indices[j]]);

					if (verticeSegment[indices[j]]!=currentSegmentID)
						LOGERROR(wxT("<GeometryTriangles>Divide: Wrong segment"));
					((GeometryTriangles*)segments[currentSegmentID])->indices.push_back(verticeMapping[indices[j]]);
				}
			}
			LOGDEBUG(wxString::Format(wxT("Divide %lu"), (unsigned long)segments.size()));
			return;
		}

		///Divide geometry to primitives, geometry is empty or has only one primitive then it returns empty list
		virtual void DivideToPrimitives(std::vector<Geometry*> &segments) const
		{
			segments.clear();
			if (vertices.empty())
				return;
			std::vector<size_t> verticeSegment(vertices.size(), (size_t)-1);
			std::vector<std::vector<size_t> > segmentVertices(indices.size());
			size_t currentSegmentID = 0, maxSegmentID = 0;
			for (size_t i = 0; i < triDescList.size(); ++i)
			{
				size_t startIndex = triDescList[i].first;
				size_t endIndex = (i + 1 < triDescList.size()) ? triDescList[i + 1].first : indices.size();
				size_t step = endIndex - startIndex;
				switch (triDescList[i].second)
				{
				case GeometryTriangles::TRI_LIST:
					step = 3;
					break;
				case GeometryTriangles::LINES:
					step = 2;
					break;
				case GeometryTriangles::POINTS:
					step = 1;
					break;
				case GeometryTriangles::QUADS:
					step = 4;
					break;
				default:
					break;
				}
				segments.push_back(new GeometryTriangles(GEOMETRY_ID_NONE));
				((GeometryTriangles*)segments.back())->triDescList.push_back(std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE>(0, triDescList[i].second));
				for (size_t j = startIndex; j < endIndex; ++j)
				{
					if (j - startIndex >= step)
					{
						segments.push_back(new GeometryTriangles(GEOMETRY_ID_NONE));
						((GeometryTriangles*)segments.back())->triDescList.push_back(std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE>(0, triDescList[i].second));
						startIndex = j;
					}
					((GeometryTriangles*)segments.back())->vertices.push_back(vertices[indices[j]]);
					if (normals.size() > indices[j])
						((GeometryTriangles*)segments.back())->normals.push_back(normals[indices[j]]);
					if (texCoords.size() > indices[j])
						((GeometryTriangles*)segments.back())->texCoords.push_back(texCoords[indices[j]]);
					((GeometryTriangles*)segments.back())->indices.push_back(j - startIndex);
				}
			}
			if (segments.size() == 1)
			{
				delete segments[0];
				segments.clear();
			}
			LOGDEBUG(wxString::Format(wxT("DivideToPrimitives %lu"), (unsigned long)segments.size()));
			return;
		}

		virtual void FlipNormals()
		{
			for (auto &normal : normals)
				normal=-normal;
		}

		virtual void TextureMappingPlanar(const Matrix &texMatrix)
		{
			texCoords.resize(vertices.size());
			for (size_t i = 0; i < vertices.size(); ++i)
			{
				Vector3 point(texMatrix.TransformCoord(vertices[i]));
				texCoords[i]=Vector2(point.x, point.y);
			}
		}

		//it is not fully general solution, moreover lines and points will be lost
		virtual void TextureMappingTriplanar(const Matrix &texMatrix)
		{
			std::vector<Index> inds;
			Triangulate(inds);
			texCoords.clear();
			texCoords.reserve(inds.size()*3);
			std::vector<Vector3> verticesTmp;
			verticesTmp.reserve(inds.size()*3);
			std::vector<Vector3> normalsTmp;
			normalsTmp.reserve(inds.size()*3);
			indices.clear();
			indices.reserve(inds.size());
			triDescList.clear();
			std::pair<GeometryTriangles::Index, GeometryTriangles::PRIMITIVE_TYPE> indType(0, GeometryTriangles::TRI_LIST);
			triDescList.push_back(indType);
			for (size_t i = 0; i < inds.size(); i+=3)
			{
				Vector3 &v1=vertices[inds[i]];
				Vector3 &v2=vertices[inds[i+1]];
				Vector3 &v3=vertices[inds[i+2]];
				verticesTmp.push_back(v1);
				verticesTmp.push_back(v2);
				verticesTmp.push_back(v3);
				Vector3 normalTr(1,0,0);
				if ((v1!=v2) && (v2!=v3) && (v3!=v1))
					normalTr=((v2-v1).Cross(v3-v1)).Normalize();
				//normalTr=texMatrix.TransformNormal(normalTr);/unnecessary, the transformation and comparison work without it
				normalTr[0]=fabs(normalTr[0]);
				normalTr[1]=fabs(normalTr[1]);
				normalTr[2]=fabs(normalTr[2]);
				Vector3 point1(texMatrix.TransformCoord(v1));
				Vector3 point2(texMatrix.TransformCoord(v2));
				Vector3 point3(texMatrix.TransformCoord(v3));
				if ((normalTr[0]>normalTr[1]) && (normalTr[0]>normalTr[2])) //X plane
				{
					texCoords.push_back(Vector2(point1.y, point1.z));
					texCoords.push_back(Vector2(point2.y, point2.z));
					texCoords.push_back(Vector2(point3.y, point3.z));
				}
				else if (normalTr[1]>normalTr[2]) //Y plane
				{
					texCoords.push_back(Vector2(point1.x, point1.z));
					texCoords.push_back(Vector2(point2.x, point2.z));
					texCoords.push_back(Vector2(point3.x, point3.z));
				}
				else //Z plane
				{
					texCoords.push_back(Vector2(point1.x, point1.y));
					texCoords.push_back(Vector2(point2.x, point2.y));
					texCoords.push_back(Vector2(point3.x, point3.y));
				}
				if (normals.size()>inds[i])
				{
					normalsTmp.push_back(normals[inds[i]]);
					normalsTmp.push_back(normals[inds[i+1]]);
					normalsTmp.push_back(normals[inds[i+2]]);
				}
				indices.push_back(i);
				indices.push_back(i+1);
				indices.push_back(i+2);
			}
			vertices.assign(verticesTmp.begin(), verticesTmp.end());
			normals.assign(normalsTmp.begin(), normalsTmp.end());
		}

		virtual void TextureMappingSpherical(const Matrix &texMatrix)
		{
			texCoords.resize(vertices.size());
			Matrix rotation(texMatrix.ExtractRotation());
			Vector3 position(texMatrix.ExtractTranslation());
			for (size_t i = 0; i < vertices.size(); ++i)
			{
				Vector3 point(rotation.TransformCoord((vertices[i]-position).Normalize()));
				texCoords[i]=Vector2((atan2(point.y, -point.x)+M_PI)/(2.0f*M_PI), acosf(-point.z)/M_PI);
			}
		}

		virtual void Clear()
		{
			vertices.clear();
			indices.clear();
			normals.clear();
			texCoords.clear();
			triDescList.clear();
			nurbsGeometry.clear();
		}

		struct VerticesVectorSort2
		{
			std::vector<Vector3> *verts;
			std::vector<Vector3> *norms;
			std::vector<GeometryTriangles::TexCoord> *texs;
			VerticesVectorSort2(std::vector<Vector3> *_verts, std::vector<Vector3> *_norms, std::vector<GeometryTriangles::TexCoord> *_texs): verts(_verts), norms(_norms), texs(_texs) {}
			bool operator()(size_t rpStart, size_t rpEnd)
			{
				if (verts->at(rpStart) < verts->at(rpEnd))
					return true;
				else if (verts->at(rpStart) == verts->at(rpEnd))
				{
					if (norms->size() == verts->size())
					{
						if (norms->at(rpStart) < norms->at(rpEnd))
							return true;
						else if (norms->at(rpStart) == norms->at(rpEnd))
						{
							if (texs->size() == verts->size())
								if (texs->at(rpStart) < texs->at(rpEnd))
									return true;
						}
					}
					else if (texs->size() == verts->size())
					{
						if (texs->at(rpStart) < texs->at(rpEnd))
							return true;
					}
				}
				return false;
			}
		};

		size_t UnifyVertices()
		{
			//sort vertices
			std::vector<size_t> ind(vertices.size());
			for (size_t i=0; i<ind.size(); i++)
				ind[i]=GeometryTriangles::Index(i);
			//std::sort(ind.begin(), ind.end(), VerticesVectorSort(vertices));
			std::sort(ind.begin(), ind.end(), VerticesVectorSort2(&vertices, &normals, &texCoords));
			//test equivalence and merge
			bool geomChanged = false;
			std::vector<std::vector<size_t> > vertsref(vertices.size());
			for (size_t i=0; i<indices.size(); i++)
				vertsref[indices[i]].push_back(i);
			size_t lastNotMerged = 0;
			for (size_t i=1; i<ind.size(); i++)
			{
				if (vertices[ind[lastNotMerged]]!=vertices[ind[i]])
				{
					lastNotMerged = i;
					continue;
				}
				if ((normals.size()==vertices.size()) && (normals[ind[lastNotMerged]]!=normals[ind[i]]))
				{
					lastNotMerged = i;
					continue;
				}
				if ((texCoords.size()==vertices.size()) && (texCoords[ind[lastNotMerged]]!=texCoords[ind[i]]))
				{
					lastNotMerged = i;
					continue;
				}
				//replace in indices
				while (!vertsref[ind[i]].empty())
				{
					indices[vertsref[ind[i]].back()]=(GeometryTriangles::Index)ind[lastNotMerged];
					vertsref[ind[lastNotMerged]].push_back(vertsref[ind[i]].back());
					vertsref[ind[i]].pop_back();
				}
				geomChanged = true;
			}
			if (!geomChanged)
			{
				for (auto &verts : vertsref)
				{
					if (verts.empty())
					{
						geomChanged=true;
						break;
					}
				}
			}
			if (geomChanged)
			{
				size_t uniCount=0;
				//remove unused vertices
				size_t lastUsed = vertsref.size()-1;
				for (size_t i=0; i<vertsref.size() && i<=lastUsed; i++)
				{
					if (vertsref[i].empty())
					{
						while ((lastUsed<vertsref.size()) && vertsref[lastUsed].empty())
						{
							lastUsed--;
							if (normals.size()==vertices.size())
								normals.pop_back();
							if (texCoords.size()==vertices.size())
								texCoords.pop_back();
							vertices.pop_back();
							uniCount++;
						}
						if ((i>=lastUsed) || (i>=vertices.size()))
							break;
						vertices[i]=vertices[lastUsed];
						if (normals.size()==vertices.size())
							normals[i]=normals[lastUsed];
						if (texCoords.size()==vertices.size())
							texCoords[i]=texCoords[lastUsed];
						while (vertsref[lastUsed].size())
						{
							indices[vertsref[lastUsed].back()]=(GeometryTriangles::Index)i;
							vertsref[i].push_back(vertsref[lastUsed].back());
							vertsref[lastUsed].pop_back();
						}
					}
				}
				return uniCount;
			}
			return 0;
		}

		///Get string with  properties - \c Geometry overload
		virtual wxString ToString() const
		{
			wxString res(Geometry::ToString());
			res+=wxString::Format(wxT("Vertices (%lu):"), (unsigned long)vertices.size());
			for (auto &vertex : vertices)
				res+=wxT(" ")+vertex.ToString();

			res+=wxString::Format(wxT("\nNormals (%lu):"), (unsigned long)normals.size());
			for (auto &normal : normals)
				res+=wxT(" ")+normal.ToString();

			res+=wxString::Format(wxT("\nTexture coordinates (%lu):"), (unsigned long)texCoords.size());
			for (auto &texCoord : texCoords)
				res+=wxT(" ")+texCoord.ToString();

			res+=wxString::Format(wxT("\nIndices (%lu):"), (unsigned long)indices.size());
			size_t indiceIndex = 0;
			size_t maxInd = indices.size();
			for (GeometryTriangles::TriDescList::const_iterator indPrimType = triDescList.begin(); indPrimType != triDescList.end(); ++indPrimType)
			{
				res+=wxString::Format(wxT("\n  Primitive type: %i"), indPrimType->second);
				const wxString primitiveNames[]={wxT("TRI_LIST"), wxT("TRI_FAN"), wxT("TRI_STRIP"), wxT("POLYGON"), wxT("LINES"), wxT("LINE_STRIP"), wxT("POINTS"), wxT("QUADS")};
				if (indPrimType->second<=(int)QUADS)
					res+=wxT(" (")+primitiveNames[indPrimType->second]+wxT(")\n");
				else
					res+=wxT(" (unknown - fix me)\n");
				GeometryTriangles::TriDescList::const_iterator endIndPrimType = indPrimType + 1;
				size_t endInd = ( endIndPrimType == triDescList.end() ? maxInd : endIndPrimType->first );
				endInd = __min( endInd, maxInd );

				for ( ;	indiceIndex < endInd; indiceIndex++)
					res+=wxString::Format(wxT(" %i"), (unsigned int)indices[indiceIndex]);
			}
			res+=wxT("\n");
			res+=wxString::Format(wxT("NURBS: count=%lu\n"), (unsigned long)nurbsGeometry.size());
			for (size_t i=0; i<nurbsGeometry.size(); i++)
				res+=wxString::Format(wxT("%lu. "), (unsigned long)i)+nurbsGeometry[i].ToString()+wxT("\n");
			return res;
		}

		///Get copy of instance - \c Geometry overload
		virtual Geometry * Clone() const
		{
			return new GeometryTriangles(*this);
		}

		private:
			GeometryTriangles &operator =(const GeometryTriangles & gt);
	};
}

#endif
