#include <vector>
#include <iostream>
#include <algorithm>
#include <string>
#include <fstream>
#include <stack>
#include <set>
#include <map>

#include <unistd.h>

#include <glm/vec3.hpp>
#include <glm/gtc/epsilon.hpp>
#include <glm/gtc/random.hpp>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
#include <OpenMesh/Tools/Decimater/DecimaterT.hh>
#include <OpenMesh/Tools/Decimater/ModProgMeshT.hh>
#include <OpenMesh/Tools/Decimater/ModQuadricT.hh>
#include <OpenMesh/Tools/Decimater/ModAspectRatioT.hh>
#include <OpenMesh/Tools/Decimater/ModEdgeLengthT.hh>

#include "PerlinNoise.hpp"
#include "cola.hpp"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"
#include "RTree.h"
// #include "MarchingCube.h"
#include "mersenne-twister.h"
#include "LSCM.h"

using namespace std;
using namespace glm;
using namespace siv;
using namespace OpenMesh;

struct MyTraits: public OpenMesh::DefaultTraits {
	VertexAttributes( OpenMesh::Attributes::Normal |
			OpenMesh::Attributes::TexCoord2D );
};

typedef TriMesh_ArrayKernelT<MyTraits> TriMesh;
typedef Decimater::DecimaterT<TriMesh> TriDecimater;
typedef Decimater::ModProgMeshT<TriMesh>::Handle HModProgMesh;
typedef Decimater::ModQuadricT<TriMesh>::Handle HModQuadric;
typedef Decimater::ModAspectRatioT<TriMesh>::Handle HModAspectRatio;
typedef Decimater::ModEdgeLengthT<TriMesh>::Handle HModEdgeLength;

class progress_callback
{
	void begin(vector<string> steps);
	void beginStep(int iStep, vector<string> substeps);
	void progress(int iSubStep, float percent);
	void endStep();
	void end();
};

class vertices {
public:
	vector<vec3> pos;
	vector<vec3> norm;
	vector<vec3> uv;

	void reserve(size_t n) {
		pos.reserve(n);
		norm.reserve(n);
		uv.reserve(n);
	}

	int append(const vertices &m) {
		size_t idx = pos.size();
		pos.insert(pos.end(), m.pos.begin(), m.pos.end());
		norm.insert(norm.end(), m.norm.begin(), m.norm.end());
		uv.insert(uv.end(), m.uv.begin(), m.uv.end());
		return idx;
	}

	int append(const vec3 &p, const vec3 &n, const vec3 &u) {
		size_t idx = pos.size();
		pos.push_back(p);
		norm.push_back(n);
		uv.push_back(u);
		return idx;
	}

};

struct edge {
	int v[2];
	int t[2];

	edge(int v1, int v2, int t1) {
		v[0] = v1;
		v[1] = v2;
		if (v1 > v2)
			swap(v[0], v[1]);

		t[0] = t1;
		t[1] = -1;
	}

	edge(int v1, int v2) {
		v[0] = v1;
		v[1] = v2;

		if (v1 > v2)
			swap(v[0], v[1]);

		t[0] = -1;
		t[1] = -1;
	}
};

struct triangle {
	int v[3];
	//int e[3];
};

bool operator<(const triangle &lhs, const triangle &rhs) {
	if (lhs.v[0] < rhs.v[0])
		return true;
	else if (lhs.v[0] > rhs.v[0])
		return false;
	else {
		if (lhs.v[1] < rhs.v[1])
			return true;
		else if (lhs.v[1] > rhs.v[1])
			return false;
		else {
			if (lhs.v[2] < rhs.v[2])
				return true;
			else if (lhs.v[2] > rhs.v[2])
				return false;
			else {
				return false;
			}

		}
	}
}

struct quant {
//	int v[3];
//
//	quant(float a, float b, float c) {
//		v[0] = 1e4 * a;
//		v[1] = 1e4 * b;
//		v[2] = 1e4 * c;
//	}

	float v[3];

	quant(float a, float b, float c) {
		v[0] = a;
		v[1] = b;
		v[2] = c;
	}
};

bool operator<(const quant &lhs, const quant &rhs) {
	if (lhs.v[0] < rhs.v[0])
		return true;
	else if (lhs.v[0] > rhs.v[0])
		return false;
	else {
		if (lhs.v[1] < rhs.v[1])
			return true;
		else if (lhs.v[1] > rhs.v[1])
			return false;
		else {
			if (lhs.v[2] < rhs.v[2])
				return true;
			else if (lhs.v[2] > rhs.v[2])
				return false;
			else {
				return false;
			}

		}
	}
}

bool operator<(const edge &lhs, const edge &rhs) {
	if (lhs.v[0] < rhs.v[0])
		return true;
	else if (lhs.v[0] > rhs.v[0])
		return false;
	else {
		if (lhs.v[1] < rhs.v[1])
			return true;
		else if (lhs.v[1] > rhs.v[1])
			return false;
		else {
			return false;
		}
	}
}

void save(ostream &out, const vertices &verts) {
	for (size_t i = 0; i < verts.pos.size(); i++) {
		out << "v " << verts.pos[i].x << " " << verts.pos[i].y << " "
				<< verts.pos[i].z << "\n";
	}
	for (size_t i = 0; i < verts.norm.size(); i++) {
		out << "vn " << verts.norm[i].x << " " << verts.norm[i].y << " "
				<< verts.norm[i].z << "\n";
	}
	for (size_t i = 0; i < verts.uv.size(); i++) {
		out << "vt " << verts.uv[i].x << " " << verts.uv[i].y << " "
				<< verts.uv[i].z << "\n";
	}
}

void save_ascii(ostream &out, const vertices &verts) {
	for (size_t i = 0; i < verts.pos.size(); i++) {
		out << verts.pos[i].x << " " << verts.pos[i].y << " " << verts.pos[i].z
				<< " ";
		out << verts.norm[i].x << " " << verts.norm[i].y << " "
				<< verts.norm[i].z << "\n";
	}
}

void save(ostream &out, const vector<triangle> &tris) {
	for (size_t i = 0; i < tris.size(); i++) {
		int v0 = tris[i].v[0] + 1;
		int v1 = tris[i].v[1] + 1;
		int v2 = tris[i].v[2] + 1;
		out << "f " << v0 << "/" << v0 << "/" << v0 << " ";
		out << v1 << "/" << v1 << "/" << v1 << " ";
		out << v2 << "/" << v2 << "/" << v2 << "\n";
	}
}

struct closest_vertex_context {
	int ignore[3];
	int v;
	float d;
	vec3 p;
	const vector<vec3> &pos;
	closest_vertex_context(const vector<vec3> &pos, const vec3 &p) :
			p(p), pos(pos) {
		v = -1;
		d = 1;
	}
};

// find closest vertex not part of first triangle
bool closest_vertex(int v, void *data) {
	closest_vertex_context *context = (closest_vertex_context *) data;

	for (size_t i = 0; i < 3; i++) {
		if (context->ignore[i] == v) {
			return true;
		}
	}

	// found closest vertex, add edges and triangle
//	triangle tr;
//	tr.v[0] = context->v0;
//	tr.v[1] = context->v1;

	float d = distance(context->p, context->pos[v]);
	if (d < context->d) {
//		tr.v[2] = v;
//		if (context->triset.find(tr) == context->triset.end()) {
		context->d = d;
		context->v = v;
//		}
	}
	return true;
}

// http://blackpawn.com/texts/pointinpoly/
bool SameSide(const vec3 &p1, const vec3 &p2, const vec3 &a, const vec3 &b) {
	vec2 cp1 = cross(b - a, p1 - a);
	vec2 cp2 = cross(b - a, p2 - a);
	if (dot(cp1, cp2) >= 0)
		return true;
	else
		return false;
}

bool PointInTriangle(const vec3 &p, const vec3 &a, const vec3 &b,
		const vec3 &c) {
	if (SameSide(p, a, b, c) && SameSide(p, b, a, c) && SameSide(p, c, a, b))
		return true;
	else
		return false;
}

int Main(int argc, char* argv[]) {
	cola::parser parser;
	parser.define("help", "Show this help").alias("h");
	parser.define("reduce", "Reduce mesh").alias("r").with_arg<int>(50000);
	parser.define("control", "Print control image").alias("c");
	parser.define("frequency", "Noise frequency").alias("f").with_arg<cola::real>(
			2);
	parser.define("octave", "Noise octaves").alias("o").with_arg<int>(2);
	parser.define("meteroids", "Meteroid craters").alias("m").with_arg<int>(
			1000);
	parser.define("seed", "Seed for noise").alias("s").with_arg<int>(time(0));
	parser.define("threshold", "Threshold for borders").alias("t").with_arg<
			cola::real>(0.2);
	parser.define("mingrad", "Minimum gradient").alias("g").with_arg<cola::real>(
			4);
	parser.define("pixels", "Number of pixels").alias("p").with_arg<int>(128);
	parser.define("texture-size", "Number of pixels").with_arg<int>(4096);
	parser.define("save-points", "Save points");
	parser.define("scale-x", "Scale in x").with_arg<cola::real>(-1);
	parser.define("scale-y", "Scale in y").with_arg<cola::real>(-1);
	parser.define("scale-z", "Scale in z").with_arg<cola::real>(-1);

	parser.parse(argc, argv);

	if (parser.is_passed("help")) {
		parser.easy_usage("generate random astroids");
		return 0;
	}

	bool show_progress = !isatty(fileno(stdin));

// 3d noise, exp falloff
	const size_t n = parser.get<int>("pixels");
	const double frequency = parser.get<cola::real>("frequency");
	const int octave = parser.get<int>("octave");
	const double t = parser.get<cola::real>("threshold");
	const double mingrad = parser.get<cola::real>("mingrad");
	const int seed = parser.get<int>("seed");
	const int meteroids = parser.get<int>("meteroids");
	const int texture_size = parser.get<int>("texture-size");

	srand(seed);

	const float scaleX =
			parser.is_passed("scale-x") ?
					parser.get<cola::real>("scale-x") : 1. + randf_cc() * 4;
	const float scaleY =
			parser.is_passed("scale-y") ?
					parser.get<cola::real>("scale-y") : 1. + randf_cc() * 4;
	const float scaleZ =
			parser.is_passed("scale-z") ?
					parser.get<cola::real>("scale-z") : 1. + randf_cc() * 4;
	const bool reduce = parser.is_passed("reduce");
	const int reduceTo = parser.get<int>("reduce");

	cout << "Generate noise field" << endl;
	cout << " frequency: " << frequency << endl;
	cout << " octave:    " << octave << endl;
	cout << " threshold: " << t << endl;
	cout << " pixels:    " << n << endl;
	cout << " seed:      " << seed << endl;
	cout << " scale:     " << scaleX << " " << scaleY << " " << scaleZ << endl;

	const size_t nx = n * n;
	const size_t ny = n;
	const size_t nz = 1;

	PerlinNoise noise(seed);
	const double s = (1.0d / n);
	const double o = 0.5;
	const double minScale = std::min(scaleZ, std::min(scaleX, scaleY));
	vector<float> voxels(n * n * n);
	for (size_t ix = 0; ix < n; ix++) {
		if (show_progress)
			(cout << "\r  " << int(100. * float(ix) / n) << "%").flush();
		double x = ix * s - o;
		size_t ox = ix * n * n;
		for (size_t iy = 0; iy < n; iy++) {
			double y = iy * s - o;
			size_t oy = ox + iy * n;
			for (size_t iz = 0; iz < n; iz++) {
				double z = iz * s - o;
				float v = (float) noise.octaveNoise0_1(x * frequency,
						y * frequency, z * frequency, octave);
				v *= std::max(
						1.
								- 4. / minScale
										* (scaleX * x * x + scaleY * y * y
												+ scaleZ * z * z), 0.);
//				if (v > t)
//					v = 1;
//				else
//					v = 0;
				voxels[oy + iz] = v;

			}
		}
	}

	// https://en.wikipedia.org/wiki/File:Vesta_Cratered_terrain_with_hills_and_ridges.jpg
	// shoot from outside, find first threshold and create hole
	for (int i = 0; i < meteroids; i++) {
		vec3 dir = sphericalRand(1.0f);
		vec3 pos = vec3(-0.7) * dir;
		vec3 idx = (pos + vec3(0.5)) * vec3(n);

		bool hit = false;
		for (size_t j = 0; !hit && j < (2 * n); j++) {
			pos += dir * vec3(s);
			idx = (pos + vec3(0.5)) * vec3(n);
			if (any(lessThan(idx, vec3(0)))) {
				continue;
			}
			if (any(greaterThan(idx, vec3(n)))) {
				continue;
			}

			float v =
					voxels[nx * int(idx.x) + ny * int(idx.y) + nz * int(idx.z)];
			if (v > 0.2f) {
				hit = true;
			}
		}

		if (hit) {
			float r = 0.005 + pow(randf_cc(), 4096) * 0.2;
			pos -= dir * vec3(r / 1.5f);
			idx = (pos + vec3(0.5)) * vec3(n);

			size_t nr = r / s;
			vec3 idx_min = max(vec3(0), idx - vec3(nr));
			vec3 idx_max = min(vec3(n), idx + vec3(nr));
			for (size_t ix = idx_min.x; ix < idx_max.x; ix++) {
				for (size_t iy = idx_min.y; iy < idx_max.y; iy++) {
					for (size_t iz = idx_min.z; iz < idx_max.z; iz++) {
						float d = length(vec3(ix, iy, iz) - idx) * s;
						if (d > r)
							continue;
						voxels[nx * ix + ny * iy + nz * iz] = 0;
					}
				}
			}
		} else {
			cout << "no hit!" << endl;
		}

	}

	if (parser.is_passed("control")) {
		cout << "Write control image" << endl;
		vector<uint8> img(n * n);
		for (size_t ix = 0; ix < n; ix++) {
			for (size_t iy = 0; iy < n; iy++) {
				img[ix * n + iy] = voxels[ix * nx + iy * ny + n / 2] * 255;
			}
		}

		stbi_write_bmp("control.bmp", n, n, 1, img.data());
	}

	if (parser.is_passed("save-points")) {

		cout << "Create points from voxels" << endl;
		cout << "  Minimum gradient: " << mingrad << endl;

		vec3 vecs[3 * 3 * 3];
		for (int ix = 0; ix < 3; ix++) {
			for (int iy = 0; iy < 3; iy++) {
				for (int iz = 0; iz < 3; iz++) {
					vecs[ix * 9 + iy * 3 + iz] = normalize(
							vec3(ix - 1, iy - 1, iz - 1));
				}
			}
		}

// center to zero
		vecs[9 + 3 + 1] = vec3(0, 0, 0);

		const size_t nx = n * n;
		const size_t ny = n;
		const size_t nz = 1;

// create vertices with normal for each border pixel
		vertices verts;

		for (size_t ix = 1; ix < (n - 1); ix++) {
			double x = ix * s - o;
			size_t ox = ix * n * n;
			for (size_t iy = 1; iy < (n - 1); iy++) {
				double y = iy * s - o;
				size_t oy = ox + iy * n;
				for (size_t iz = 1; iz < (n - 1); iz++) {
					size_t io = oy + iz;
					if (voxels[io] == 0)
						continue;

					vec3 gradient = vec3(0, 0, 0);
					for (int jx = 0; jx < 3; jx++) {
						for (int jy = 0; jy < 3; jy++) {
							for (int jz = 0; jz < 3; jz++) {
								gradient +=
										vecs[jx * 9 + jy * 3 + jz]
												* (float) voxels[io
														+ nx * (jx - 1)
														+ ny * (jy - 1)
														+ nz * (jz - 1)];
							}
						}
					}

					// no border = no gradient
					if (length(gradient) < mingrad + 0.1)
						continue;

					double z = iz * s - o;

					verts.append(vec3(x, y, z), normalize(gradient) * -1.f,
							vec3(0, 0, 0));

				}
			}
		}

		cout << "Save points" << endl;
		ofstream ofile("points.obj");
		save(ofile, verts);

		ofstream oafile("points.ascii");
		save_ascii(oafile, verts);
	}

	cout << "Marching cubes" << endl;

	TriMesh mesh;

	vertices verts;
	vector<triangle> tris;
	map<quant, int> vertexLookup;
	RTree<int, float, 3, float> tree;

	size_t nT = 0;
	for (size_t ix = 0; ix < (n - 1); ix++) {
		if (show_progress)
			(cout << "\r  " << int(100. * float(ix) / n) << "%").flush();

		double x = ix * s - o;
		size_t ox = ix * n * n;
		for (size_t iy = 0; iy < (n - 1); iy++) {
			double y = iy * s - o;
			size_t oy = ox + iy * n;
			for (size_t iz = 0; iz < (n - 1); iz++) {
				size_t io = oy + iz;
				double z = iz * s - o;
//
//				GRIDCELL gridCell;
//				gridCell.p[0] = vec3(x, y, z);
//				gridCell.p[1] = vec3(x + s, y, z);
//				gridCell.p[2] = vec3(x + s, y + s, z);
//				gridCell.p[3] = vec3(x, y + s, z);
//				gridCell.p[4] = vec3(x, y, z + s);
//				gridCell.p[5] = vec3(x + s, y, z + s);
//				gridCell.p[6] = vec3(x + s, y + s, z + s);
//				gridCell.p[7] = vec3(x, y + s, z + s);
//
//				gridCell.val[0] = voxels[io];
//				gridCell.val[1] = voxels[io + nx];
//				gridCell.val[2] = voxels[io + nx + ny];
//				gridCell.val[3] = voxels[io + ny];
//				gridCell.val[4] = voxels[io + nz];
//				gridCell.val[5] = voxels[io + nx + nz];
//				gridCell.val[6] = voxels[io + nx + ny + nz];
//				gridCell.val[7] = voxels[io + ny + nz];
//				TRIANGLE triangles[6];
//				int nTris = Polygonise(gridCell, t, triangles);
//				nT += nTris;
//				for (int iTri = 0; iTri < nTris; iTri++) {
//					vec3 *ps = triangles[iTri].p;
//					if (all(epsilonEqual(ps[0], ps[1], 1e-8f))
//							|| all(epsilonEqual(ps[0], ps[2], 1e-8f))
//							|| all(epsilonEqual(ps[2], ps[1], 1e-8f)))
//						continue;
//					TriMesh::VertexHandle vhandle[3];
//					for (size_t k = 0; k < 3; k++) {
//						vec3 &p = ps[k];
//						quant q(p.x, p.y, p.z);
//						if (reduce) {
//							map<quant, int>::iterator vit = vertexLookup.find(
//									q);
//							if (vit != vertexLookup.end()) {
//								TriMesh::VertexHandle h = mesh.vertex_handle(
//										vit->second);
//								vhandle[k] = h;
//							}
//						}
//						if (!vhandle[k].is_valid()) {
//							TriMesh::Point pnt(p.x, p.y, p.z);
//							vhandle[k] = mesh.add_vertex(pnt);
//							if (reduce)
//								vertexLookup[q] = vhandle[k].idx();
//
//						}
//					}
//
//					mesh.add_face(vhandle, 3);
//				}
			}
			//(cout << ".").flush();
		}
		//(cout << " " << nT << " " << vertexLookup.size() << "\n").flush();
	}
#if 0

// create triangles
	RTree<int, float, 3, float> tree;
	for (size_t i = 0; i < verts.pos.size(); i++) {
		float p[3] = {verts.pos[i].x, verts.pos[i].y, verts.pos[i].z};
		tree.Insert(p, p, i);
	}

// find closest vertex to first
	int second_vertex = -1;
	float second_dist = 1.0;
	for (size_t j = 1; j < verts.pos.size(); j++) {
		float d = distance(verts.pos[j], verts.pos[0]);
		if (d < second_dist) {
			second_dist = d;
			second_vertex = j;
		}
	}

	cout << "Create triangles from vertices" << endl;
	vector<triangle> tris;
	vector<edge> edges;
	stack<int> uedges;

	edges.push_back(edge(0, second_vertex));
	uedges.push(edges.size() - 1);

	while (uedges.size() > 0) {
		int ei = uedges.top();
		int v0 = edges[ei].v[0];
		int v1 = edges[ei].v[1];
		int t0 = edges[ei].t[0];
		int t1 = edges[ei].t[1];

		vec3 center = (verts.pos[v0] + verts.pos[v1]) / 2.f;
		closest_vertex_context context(verts.pos, center);

// first triangle set?
		if (t0 >= 0) {
			// remove from unfinished
			uedges.pop();

			// first triangle
			for (size_t i = 0; i < 3; i++)
			context.ignore[i] = tris[t0].v[i];
		} else {
			// change first triangle
			context.ignore[0] = v0;
			context.ignore[1] = v1;
			context.ignore[2] = v1;
		}

// find closest vertex not in edge of first triangle
		float range = second_dist;
		while (context.v < 0) {
			range *= 2;
			float smin[3] = {center.x - range, center.y - range, center.z
				- range};
			float smax[3] = {center.x + range, center.y + range, center.z
				+ range};
			if (tree.Search(smin, smax, closest_vertex, &context) == 0)
			cout << "no results!" << endl;
		}

// found closest vertex, add edges and triangle
		triangle tr;
		tr.v[0] = v0;
		tr.v[1] = v1;
		tr.v[2] = context.v;

		int e1 = edges.size();
		int e2 = e1 + 1;
		int t = tris.size();

		if (t0 < 0) {
			edges[ei].t[0] = t;
		} else {
			edges[ei].t[1] = t;
		}

		edges.push_back(edge(v0, context.v, t));
		edges.push_back(edge(v1, context.v, t));
		uedges.push(e1);
		uedges.push(e2);

//		tr.e[0] = ei;
//		tr.e[1] = e1;
//		tr.e[2] = e2;
		tris.push_back(tr);
		triset.insert(tr);

		if (tris.size() % 1000 == 0) {
			cout << tris.size() << " ";
			cout << uedges.size() << endl;
		}
	}

#endif

	cout << "Stage1 Vertices: " << mesh.n_vertices() << endl;
	cout << "Stage1 Faces: " << mesh.n_faces() << endl;
	cout << "Stage1 Save file" << endl;
	OpenMesh::IO::write_mesh(mesh, "stage1.obj");

	cout << "Stage 2: Smooth mesh" << endl;

	{
// this vector stores the computed centers of gravity
		vector<TriMesh::Point> cogs;
		vector<TriMesh::Point>::iterator cog_it;
		cogs.reserve(mesh.n_vertices());
// smoothing mesh argv[1] times
		TriMesh::VertexIter v_it, v_end(mesh.vertices_end());
		TriMesh::VertexVertexIter vv_it;
		TriMesh::Point cog;
		TriMesh::Scalar valence;
		for (int i = 0; i < 2; ++i) {
			cogs.clear();
			for (v_it = mesh.vertices_begin(); v_it != v_end; ++v_it) {
				cog[0] = cog[1] = cog[2] = valence = 0.0;

				for (vv_it = mesh.vv_iter(*v_it); vv_it.is_valid(); ++vv_it) {
					cog += mesh.point(*vv_it);
					++valence;
				}
				cogs.push_back(cog / valence);
			}

			for (v_it = mesh.vertices_begin(), cog_it = cogs.begin();
					v_it != v_end; ++v_it, ++cog_it)
				if (!mesh.is_boundary(*v_it))
					mesh.set_point(*v_it, *cog_it);
		}

	}

	cout << "Stage2: Vertices: " << mesh.n_vertices() << endl;
	cout << "Stage2: Faces: " << mesh.n_faces() << endl;
	cout << "Stage2: Save file" << endl;
	OpenMesh::IO::write_mesh(mesh, "stage2.obj");

	cout << "Stage3: Optimize mesh" << endl;
#if 1
	mesh.request_face_status();
	mesh.request_edge_status();
	mesh.request_vertex_status();

	TriDecimater decimater(mesh);  // a decimater object, connected to a mesh

#if 1
	HModQuadric hModQuadric;      // use a quadric module
	decimater.add(hModQuadric);      // register module at the decimater
	cout << decimater.module(hModQuadric).name() << endl;      // module access
	/*
	 * since we need exactly one priority module (non-binary)
	 * we have to call set_binary(false) for our priority module
	 * in the case of HModQuadric, unset_max_err() calls set_binary(false) internally
	 */
	decimater.module(hModQuadric).unset_max_err();

#endif
#if 0

	HModAspectRatio hModAspectRatio;      // use a quadric module
	decimater.add(hModAspectRatio);// register module at the decimater
//decimater.module(hModAspectRatio).set_aspect_ratio(5);
	decimater.module(hModAspectRatio).set_binary(false);
	cout << decimater.module(hModAspectRatio).name() << endl;// module access

#endif
#if 0

	HModEdgeLength hModEdgeLength;      // use a quadric module
	decimater.add(hModEdgeLength);// register module at the decimater
//decimater.module(hModAspectRatio).set_aspect_ratio(5);
	decimater.module(hModEdgeLength).set_binary(false);
	cout << decimater.module(hModEdgeLength).name() << endl;// module access

#endif

//	HModProgMesh hModProgMesh;      // use a quadric module
//	decimater.add(hModProgMesh); // register module at the decimater
//	cout << decimater.module(hModProgMesh).name() << endl; // module access

	if (!decimater.initialize())
		cerr << "Init failed!" << endl;
	cout << "Stage3: Decimated " << decimater.decimate_to(reduceTo) << endl;
// after decimation: remove decimated elements from the mesh
	mesh.garbage_collection();

	cout << "Stage3: Vertices: " << mesh.n_vertices() << endl;
	cout << "Stage3: Faces: " << mesh.n_faces() << endl;
	cout << "Stage3: Save file" << endl;
	OpenMesh::IO::write_mesh(mesh, "stage3.obj");

#endif // OPT

	cout << "Stage4: UV map" << endl;
	TriMesh mappedMesh;
	float texf = 1. / texture_size;
	vector<uint8[3]> diff(texture_size * texture_size * 3);
	mesh.request_vertex_normals();
	mesh.request_face_normals();
	mesh.update_face_normals();
	mesh.update_vertex_normals();
	int rects = ceil(sqrt(mesh.n_faces())), rx = 0, ry = 0;
	float rectSize = floor(float(texture_size) / rects) / texture_size;
	size_t nfaces = mesh.n_faces(), iface = 0;
	for (auto i = mesh.faces_begin(); i != mesh.faces_end(); i++) {
		if (show_progress)
			(cout << "\r  " << int(100. * float(iface++) / nfaces) << "%").flush();

		vec3 p[3], n[3];
		int k = 0;
		TriMesh::VertexHandle h[3];
		for (auto j = mesh.fv_begin(*i); j != mesh.fv_end(*i); j++) {
			if (!j->is_valid())
				cout << "invalid vertex" << endl;
			TriMesh::Point pt = mesh.point(*j);
			p[k] = vec3(pt[0], pt[1], pt[2]);
			TriMesh::Normal nm = mesh.calc_vertex_normal(*j);
			n[k] = vec3(nm[0], nm[1], nm[2]);
			h[k] = mappedMesh.add_vertex(pt);
			mappedMesh.set_normal(h[k], nm);
			k++;
			if (k > 3)
				cout << "not a triangle!" << endl;
		}
		mappedMesh.add_face(h, 3);

		vec3 ex = normalize(p[1] - p[0]);
		vec3 ey = normalize(cross(ex, p[2] - p[0]));
		vec3 ez = normalize(cross(ex, ey));
		mat3 base = mat3(ex, ey, ez);
		vec3 t[3] = { p[0] * base, p[1] * base, p[2] * base };
		vec3 lower = min(min(t[0], t[1]), t[2]);
		vec3 upper = max(max(t[0], t[1]), t[2]);
		t[0] -= lower;
		t[1] -= lower;
		t[2] -= lower;
		vec3 extent = upper - lower;
		float s = std::max(std::max(extent.x, extent.y), extent.z);
		t[0] *= 0.8 * rectSize / s;
		t[1] *= 0.8 * rectSize / s;
		t[2] *= 0.8 * rectSize / s;
		vec3 off(rx * rectSize, 0, ry * rectSize);
		t[0] += off;
		t[1] += off;
		t[2] += off;
		mappedMesh.set_texcoord2D(h[0], TriMesh::TexCoord2D(t[0].x, t[0].z));
		mappedMesh.set_texcoord2D(h[1], TriMesh::TexCoord2D(t[1].x, t[1].z));
		mappedMesh.set_texcoord2D(h[2], TriMesh::TexCoord2D(t[2].x, t[2].z));

		//if (rx == 0 && ry == 0) {
		t[0].y = 0;
		t[1].y = 0;
		t[2].y = 0;
//			cout << "ftex " << texf << endl;
//			cout << "t[0]  " << t[0].x << " " << t[0].z << endl;
//			cout << "t[1]  " << t[1].x << " " << t[1].z << endl;
//			cout << "t[2]  " << t[2].x << " " << t[2].z << endl;
//			cout << "t[0]p " << t[0].x * ntex << " " << (1 - t[0].z) * ntex << endl;
//			cout << "t[1]p " << t[1].x * ntex << " " << (1 - t[1].z) * ntex << endl;
//			cout << "t[2]p " << t[2].x * ntex << " " << (1 - t[2].z) * ntex << endl;
		// fill texture
		// fill whole rect for now
		int rpix = floor(rectSize * texture_size);

		vec3 eb = t[1] - t[0];
		float leb = length(eb);
		vec3 ec = t[2] - t[0];
		float lec = length(ec);
//		mat2 textureBase = mat2(vec2(eb.x, eb.z), vec2(ec.x, ec.z));
//		mat2 itextureBase = inverse(textureBase);
		mat3 textureBase = mat3(eb, cross(eb, ec), ec);
		mat3 itextureBase = inverse(textureBase);

//		eb = normalize(eb);
//		ec = normalize(ec);

		vec3 db = p[1] - p[0];
		vec3 dc = p[2] - p[0];

		vec3 mb = normalize(p[1] - p[0]);
		vec3 mc = normalize(p[2] - p[0]);
//		mat2 faceBase = mat2(vec2(mb.x, mb.z), vec2(mc.x, mc.z));
//		mat2 ifaceBase = inverse(faceBase);

		mat3 faceBase = mat3(db, cross(db, dc), dc);
		mat3 ifaceBase = inverse(faceBase);

		// mat2 trafo = ifaceBase * textureBase * faceBase;

		for (int ix = rx * rpix; ix < (rx + 1) * rpix; ix++) {
			for (int iy = (texture_size - (ry + 1) * rpix);
					iy < (texture_size - ry * rpix); iy++) {

//					cout << "pix " << ix << " " << iy << endl;
//					cout << "uv  " << float(ix) * texf << " " << 1.f - float(iy) * texf << endl;

				// pixel to uv space
				vec3 res = vec3(float(ix) * texf, 0.f, 1.f - float(iy) * texf)
						- t[0];
				// cout << "r   = " << r.x << " " << r.y << endl;
				res = faceBase * itextureBase * res;

				// cout << "res = " << res.x << " " << res.y << endl;
//				vec2 res = ifaceBase * textureBase * r * faceBase;
				//
//				cout << "pix " << ix << " " << res.x << " " << res.y << " " << res.z << endl;
//				float fb = dot(r, eb) / leb;
//				float fc = dot(r, ec) / lec;
				vec3 sp = p[0] + res; //res.x * mb + res.y * mc;
				float v = (float) noise.octaveNoise0_1(sp.x * frequency * 4,
						sp.y * frequency * 4, sp.z * frequency * 4, 8);
//
				v = pow(v, 3);
				// default color
				float r = 0.12 + v * 0.72, g = 0.1 + 0.7 * v, b = 0.1 + 0.7 * v;

				diff[ix + iy * texture_size][0] = r * 255;
				diff[ix + iy * texture_size][1] = g * 255;
				diff[ix + iy * texture_size][2] = b * 255;

//				diff[ix + iy * ntex] = v * 255;
//				diff[ix + iy * ntex][0] = (32 * rx) % 255;
//				diff[ix + iy * ntex][1] = (32 * ry) % 255;
//				diff[ix + iy * ntex][2] = 0; //(0.5 + sp.z) * 255;

//				float rd = 1.f; //0.7 + 0.3 * length(r) * rects;
//				diff[ix + iy * ntex][0] = rd * (0.5 + sp.x) * 255;
//				diff[ix + iy * ntex][1] = rd * (0.5 + sp.y) * 255;
//				diff[ix + iy * ntex][2] = rd * (0.5 + sp.z) * 255;

//					diff[ix + iy * ntex][0] = length(r) * rects * 255;
//					diff[ix + iy * ntex][1] = length(r) * rects * 255;
//					diff[ix + iy * ntex][2] = length(r) * rects * 255;

//									diff[ix + iy * ntex][0] = r.x * rects * 200;
//									diff[ix + iy * ntex][1] = -r.z * rects * 200;
//									diff[ix + iy * ntex][2] = 0;

//													diff[ix + iy * ntex][0] = 128 + (res.x - 0.5) * 200;
//													diff[ix + iy * ntex][1] = 128 + (res.y - 0.5) * 200;
//													diff[ix + iy * ntex][2] = 0; //res.z * 255;

			}
			//}
		}
		rx++;
		if (rx >= rects) {
			rx = 0;
			ry++;
		}

	}
//
//	for (int ix = 0; ix < ntex; ix++) {
//		for (int iy = 0; iy < ntex; iy++) {
//
//			float v = (float) noise.octaveNoise0_1(ix * texf * frequency,
//					iy * texf * frequency, 0.5 * frequency, 4);
//			diff[ix + (ntex - iy) * ntex] = 5 + v * 250;
//		}
//	}
	stbi_write_bmp("diffuse.bmp", texture_size, texture_size, 3, diff.data());

	cout << "Stage4: Vertices: " << mesh.n_vertices() << endl;
	cout << "Stage4: Faces: " << mesh.n_faces() << endl;
	cout << "Stage4: Save file" << endl;
	OpenMesh::IO::write_mesh(mappedMesh, "stage4.obj",
			IO::Options::VertexTexCoord | IO::Options::VertexNormal);
	/*
	 for (size_t ix = 0; ix < n; ix++) {
	 for (size_t iy = 0; iy < n; iy++) {
	 img[ix * n + iy] = voxels[ix * nx + iy * ny + n / 2] * 255;
	 }
	 }


	 */
#if LSCM
	nlInitialize(argc, argv);

// nlInitialize(argc, argv);
	LSCM::IndexedMesh lscmMesh;

// copy data
	const TriMesh::Point *points = mesh.points();
	for (size_t i = 0; i < mesh.n_vertices(); i++) {
		lscmMesh.add_vertex(
				LSCM::vec3(points[i].data()[0], points[i].data()[1],
						points[i].data()[2]), LSCM::vec2(0.0, 0.0));
	}

	for (auto i = mesh.faces_begin(); i != mesh.faces_end(); i++) {
		lscmMesh.begin_facet();
		for (auto j = mesh.fv_begin(*i); j != mesh.fv_end(*i); j++) {
			lscmMesh.add_vertex_to_facet(j->idx());
		}
		lscmMesh.end_facet();
	}

	LSCM::LSCM lscm(lscmMesh);
	lscm.set_spectral(false);
	lscm.apply();

	lscmMesh.save("stage3.obj");
#endif

//	// copy uv back
//	mesh.request_vertex_texcoords2D();
//
//	if (lscmMesh.nb_vertices() != mesh.n_vertices())
//		cout << "invali" << endl;
//	for (unsigned int v = 0; v < lscmMesh.nb_vertices(); ++v) {
//		TriMesh::TexCoord2D tex(lscmMesh.vertex[v].tex_coord.x, lscmMesh.vertex[v].tex_coord.y);
//		mesh.set_texcoord2D(mesh.vertex_handle(v), tex);
//	}
//
//	cout << "Write File" << endl;
//	OpenMesh::IO::write_mesh(mesh, "stage3.obj", IO::Options::VertexTexCoord);

	return 0;
}