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

#ifndef __VCGLIB_ZONOHEDRON
#define __VCGLIB_ZONOHEDRON


#include<vcg/complex/allocate.h>
#include<map>
typedef unsigned int uint;

namespace vcg {
namespace tri {
/** \addtogroup trimesh */
//@{
    /**
				A class to build a Zonohedron.

				A zonohedron is a solid with a closed surface composed only of parallelograms.
				Given a set of input vectors (the sides of the parallelograms), it is defined
				by the convex hull of all the points which can be co

				Creates a pure-quad mesh (triangular bit-quad),
				(faces with 4 vertices are split into quads).

				USAGE:
					1) Instantiate a Zonohedron.
					2) Add input vectors at will to it, with addVector(s)
					3) When you are done, call createMesh.

    */


template <class Scalar>
class Zonohedron{
public:

	typedef Point3<Scalar> Vec3;

	Zonohedron(){}

	void addVector(Scalar x, Scalar y, Scalar z);
	void addVector(Vec3 v);
	void addVectors(const std::vector< Vec3 > );

	const std::vector< Vec3 >& vectors() const {
		return vec;
	}

	template<class MeshType>
	void createMesh( MeshType& output );

private:

	/* classes for internal use */
	/****************************/

	typedef int VecIndex; //  a number in [0..n)

	/* the signature of a vertex (a 0 or 1 per input vector) */
	struct Signature {
		std::vector< bool > v;
		Signature(){}
		Signature(int n){ v.resize(n,false); }

		bool operator == (const Signature & b) const {
			return (b.v == v);
		}
		bool operator < (const Signature & b) const {
			return (b.v < v);
		}
		Signature& set(VecIndex i, bool value){
			v[i] = value;
			return *this;
		}
		Signature& set(VecIndex i, bool valueI, VecIndex j, bool valueJ){
			v[i] = valueI;
			v[j] = valueJ;
			return *this;
		}
	};

	struct Face {
		int vert[4]; // index to vertex array
	};

	/* precomputed cross products for all pairs of vectors */
	std::vector< Vec3 > precomputedCross;

	void precompteAllCrosses(){
		precomputedCross.resize(n*n);
		for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
			precomputedCross[i*n+j] =  vec[i] ^ vec[j] ;
		}
	}

	Vec3 cross(VecIndex i, VecIndex j){
		return precomputedCross[i*n+j];
	}

	// given a vector, returns a copy pointing a unique verse
	static Vec3 uniqueVerse(Vec3 v){
		if (v.X()>0) return v;
		else if (v.X()<0) return -v;
		else if (v.Y()>0) return v;
		else if (v.Y()<0) return -v;
		else if (v.Z()>0) return v;
		return -v;
	}

	// returns signof:  (i x j) * k
	bool signOf_IxJoK(VecIndex i, VecIndex j, VecIndex k){
		bool invert = false;
		// sort i,j,k
		if (i<j) { std::swap(i,j); invert = !invert; }
		if (j<k) { std::swap(j,k); invert = !invert;
			if (i<j) { std::swap(i,j); invert = !invert; }
		}

		//Scalar res = Vec3::dot( Vec3::cross( vec[i] , vec[j] ) , vec[k] );
		Scalar res =  cross( i , j ) * vec[k] ;

		if (res==0) {
			// three coplanar vectors!
			res =  uniqueVerse( cross(i,j) ) * cross(j,k) ;
			/*if (res==0) {
				printf("COLINEAR\n");
				// three colinear vectors!
				res = Vec3::dot( uniqueVerse(vec[i]) , vec[j])*
				      Vec3::dot( uniqueVerse(vec[i]) , vec[k]);
			}*/
		}

		return ( (res>=0) != invert ); // XOR
	}

	int n; // number of input vectors
	std::vector<Vec3> vec; // input vectors

	int vertCount;
	std::vector<Face> _face;

	typedef std::map< Signature, int > VertexMap;
	VertexMap vertexMap;

	// given a vertex signature, returns index of vert (newly created or not)
	VecIndex vertexIndex(const Signature &s){
		typename VertexMap::iterator i;
		//Vec3 pos = s; //toPos(s);
		i = vertexMap.find( s );
		if (i!= vertexMap.end() ) return i->second;
		else {
			int newVertex = vertCount++;
			//vertexMap.insert(s)
			vertexMap[s] = newVertex;
			return newVertex;
		}
	}

	// given two index of vectors, returns face
	Face& face(VecIndex i, VecIndex j){
		assert(i!=j);
		assert( i*n + j < (int) _face.size() );
		return _face[i*n + j];
	}

	Vec3 toPos(const Signature &s) const{
		Vec3 res(0,0,0);
		for (int i=0; i<n; i++)
			if (s.v[i]) res += vec[i];
		return res;
	}

	void createInternalMesh() {

		n = vec.size();
		precompteAllCrosses();

		// allocate faces
		_face.resize( n*n );

		vertCount = 0;
		vertexMap.clear();

		for (int i=0; i<n; i++) {
			//::showProgress(i,n);
			for (int j=0; j<n; j++) if(i!=j)  {
				Signature s(n);
				for (int k=0; k<n; k++) if ((k!=j) && (k!=i))
				{
					s.set( k , signOf_IxJoK( i,j,k ) );
				}
				face(i,j).vert[0] = vertexIndex( s.set(i,false, j,false) );
				face(i,j).vert[1] = vertexIndex( s.set(i,false, j,true ) );
				face(i,j).vert[2] = vertexIndex( s.set(i,true,  j,true ) );
				face(i,j).vert[3] = vertexIndex( s.set(i,true,  j,false) );
			}
		}
	}


};


template<class Scalar>
void Zonohedron<Scalar>::addVectors(std::vector< Zonohedron<Scalar>::Vec3 > input){
	for (uint i=0; i<input.size(); i++) {
		addVector( input[i]);
	}
}

template<class Scalar>
void Zonohedron<Scalar>::addVector(Scalar x, Scalar y, Scalar z) {
	addVector( Vec3(x,y,z) );
}


template<class Scalar>
void Zonohedron<Scalar>::addVector(Zonohedron<Scalar>::Vec3 v){
	vec.push_back(v);
}


template<class Scalar>
template<class MeshType>
void Zonohedron<Scalar>::createMesh(MeshType &m){
	typedef MeshType Mesh;
	typedef typename Mesh::VertexPointer  MeshVertexPointer;
	typedef typename Mesh::VertexIterator MeshVertexIterator;
	typedef typename Mesh::FaceIterator   MeshFaceIterator;
	typedef typename Mesh::FaceType   MeshFace;

	createInternalMesh();

	m.Clear();
	Allocator<MeshType>::AddVertices(m,vertexMap.size());
  Allocator<MeshType>::AddFaces(m,n*(n-1) * 2);

	// assign vertex positions
	MeshVertexIterator vi=m.vert.begin();
	for (typename VertexMap::iterator i=vertexMap.begin(); i!=vertexMap.end(); i++){
		(vi + i->second )->P() = toPos( i->first );
	}

	// assegn FV connectivity
  MeshFaceIterator fi=m.face.begin();

	for (int i=0; i<n; i++) {
		for (int j=0; j<n; j++) if (i!=j) {
			const Face &f( face(i,j) );
			for (int k=0; k<2; k++) { // two tri faces per quad
				for (int w=0; w<3; w++) {
					fi->V(w) = &* (vi + f.vert[(w+k*2)%4] );
				}
				if (tri::HasPerFaceNormal(m)) {
					fi->N() = cross(i,j).normalized();
				}
				if (tri::HasPerFaceFlags(m)) {
					fi->SetF(2); // quad diagonals are faux
				}
				fi++;
			}
		}
	}


}


//@}

} // End Namespace TriMesh
} // End Namespace vcg

#endif // __VCGLIB_ZONOHEDRON