initial commit

src/AsteroidGenerator.h

@@ -0,0 +1,15 @@
+/*
+ * AsteroidGenerator.h
+ *
+ *  Created on: 24.10.2017
+ *      Author: gmueller
+ */
+
+#ifndef SRC_ASTEROIDGENERATOR_H_
+#define SRC_ASTEROIDGENERATOR_H_
+
+
+
+
+
+#endif /* SRC_ASTEROIDGENERATOR_H_ */

src/Density.h

@@ -0,0 +1,126 @@
+#pragma once
+
+// system
+#include <vector>
+
+// external
+#include <glm/vec3.hpp>
+#include <glm/gtc/random.hpp>
+
+// vendor
+#include <RTree.h>
+#include <PerlinNoise.hpp>
+#include <stb_image_write.h>
+
+struct Impact {
+	glm::vec3 pos;
+	float r;
+	bool contains(const glm::vec3 &p) const {
+		return (glm::length(pos - p) < r);
+	}
+};
+
+class Density {
+	std::vector<Impact> impacts;
+	RTree<int, float, 3, float> impactTree;
+	siv::PerlinNoise noise;
+	glm::vec3 scale;
+	float scaleFactor;
+	float frequency, octave;
+
+public:
+	Density() :
+			scale(1), scaleFactor(4), frequency(2), octave(2) {
+
+	}
+
+	void addImpact(Impact &impact) {
+		size_t idx = impacts.size();
+		impacts.push_back(impact);
+		glm::vec3 l = impact.pos - glm::vec3(impact.r);
+		glm::vec3 u = impact.pos + glm::vec3(impact.r);
+		impactTree.Insert(&l.x, &u.x, idx);
+	}
+
+	void setSeed(uint32_t seed) {
+		noise.reseed(seed);
+	}
+
+	void setScale(const glm::vec3 &scale) {
+		this->scale = scale;
+		scaleFactor = 4.f / std::min(scale.x, std::min(scale.y, scale.z));
+	}
+
+	void setFrequency(float f) {
+		frequency = f;
+	}
+
+	void setOctave(float o) {
+		octave = o;
+	}
+
+	float operator()(const glm::vec3 &p) const {
+		if (impacts.size()) {
+			bool insideImpact = false;
+			impactTree.Search(&p.x, &p.x, [&](int idx) -> bool {
+				if (impacts[idx].contains(p)) {
+					insideImpact = true;
+					return false;
+				} else {
+					return true;
+				}
+			});
+
+			if (insideImpact)
+				return 0.f;
+		}
+		float v = (float) noise.octaveNoise0_1(p.x * frequency, p.y * frequency,
+				p.z * frequency, octave);
+		glm::vec3 r = p - glm::vec3(0.5);
+		v *= std::max(1.f - scaleFactor * dot(scale * r, r), 0.f);
+		return v;
+	}
+
+	bool intersectIsolevel(const glm::vec3 &start, const glm::vec3 &end,
+			glm::vec3 &pos, float isolevel, float resolution) {
+		glm::vec3 step = glm::normalize(end - start) * glm::vec3(resolution);
+		size_t nSteps = ceil(glm::length(end - start) / resolution);
+		pos = start;
+		for (size_t iStep = 0; iStep < nSteps; iStep++) {
+			if ((*this)(pos) > isolevel)
+				return true;
+			pos += step;
+		}
+		return false;
+	}
+
+	void addRandomImpacts(size_t count, float isolevel, float resolution,
+			float minRadius, float maxRadius, int exponent) {
+		for (size_t i = 0; i < count; i++) {
+			glm::vec3 start = glm::vec3(0.5f) + glm::sphericalRand(0.6f);
+			glm::vec3 end = glm::vec3(0.5f) + glm::sphericalRand(0.6f);
+			Impact impact;
+			if (intersectIsolevel(start, end, impact.pos, isolevel,
+					resolution)) {
+				impact.r = minRadius
+						+ pow(glm::linearRand(0.f, 1.f), exponent) * maxRadius;
+				addImpact(impact);
+			}
+		}
+	}
+
+	void saveCrossSection(const char *filename, size_t pixels) {
+		float step = 1.f / pixels;
+		glm::vec3 pos(step / 2, step / 2, 0.5f);
+		std::vector<uint8_t> img(pixels * pixels);
+		for (size_t ix = 0; ix < pixels; ix++) {
+			for (size_t iy = 0; iy < pixels; iy++) {
+				glm::vec3 pos(step * (0.5f + ix), step * (0.5f + iy), 0.5f);
+				img[iy * pixels + ix] = (*this)(pos) * 255;
+			}
+		}
+
+		stbi_write_png(filename, pixels, pixels, 1, img.data(), 0);
+	}
+};
+

src/Polygoniser.h

@@ -0,0 +1,117 @@
+#pragma once
+
+#include <functional>
+#include <vector>
+#include <map>
+
+#include <glm/vec3.hpp>
+#include <glm/gtc/epsilon.hpp>
+
+#include <MarchingCube.h>
+#include <logger.h>
+
+#include "Density.h"
+
+struct CompareVec3 {
+	bool operator()(const glm::vec3 &lhs, const glm::vec3 &rhs) {
+		static const float epsilon = 1e-5;
+		if (lhs.x < rhs.x - epsilon)
+			return true;
+		else if (lhs.x > rhs.x + epsilon)
+			return false;
+		else {
+			if (lhs.y < rhs.y - epsilon)
+				return true;
+			else if (lhs.y > rhs.y + epsilon)
+				return false;
+			else {
+				if (lhs.z < rhs.z - epsilon)
+					return true;
+				else if (lhs.z > rhs.z + epsilon)
+					return false;
+				else {
+					return false;
+				}
+
+			}
+		}
+	}
+};
+
+class Polygonizer {
+public:
+	std::vector<glm::vec3> vertices;
+	std::vector<glm::uvec3> triangles;
+
+	Density &density;
+
+	Polygonizer(Density &density) :
+			density(density) {
+
+	}
+
+	void polygonize(float isolevel, float resolution) {
+		vertices.clear();
+		triangles.clear();
+
+		std::map<glm::vec3, int, CompareVec3> vertexLookup;
+
+		size_t nSteps = 1.f / resolution - 1;
+		size_t nT = 0;
+		for (size_t ix = 0; ix < nSteps; ix++) {
+			double x = (double(ix) + 0.5) * resolution;
+			for (size_t iy = 0; iy < nSteps; iy++) {
+				double y = (double(iy) + 0.5) * resolution;
+				for (size_t iz = 0; iz < nSteps; iz++) {
+					double z = (double(iz) + 0.5) * resolution;
+
+					GRIDCELL gridCell;
+					gridCell.p[0] = vec3(x, y, z);
+					gridCell.p[1] = vec3(x + resolution, y, z);
+					gridCell.p[2] = vec3(x + resolution, y + resolution, z);
+					gridCell.p[3] = vec3(x, y + resolution, z);
+					gridCell.p[4] = vec3(x, y, z + resolution);
+					gridCell.p[5] = vec3(x + resolution, y, z + resolution);
+					gridCell.p[6] = vec3(x + resolution, y + resolution,
+							z + resolution);
+					gridCell.p[7] = vec3(x, y + resolution, z + resolution);
+
+					for (size_t iCell = 0; iCell < 8; iCell++) {
+						gridCell.val[iCell] = density(gridCell.p[iCell]);
+					}
+
+					TRIANGLE tris[6];
+					int nTris = Polygonise(gridCell, isolevel, tris);
+					nT += nTris;
+					for (int iTri = 0; iTri < nTris; iTri++) {
+						vec3 *ps = tris[iTri].p;
+
+						// skip degenerate vertices
+						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;
+
+						glm::uvec3 idc;
+						for (size_t k = 0; k < 3; k++) {
+							glm::vec3 &p = tris[iTri].p[k];
+							std::map<glm::vec3, int>::iterator vit =
+									vertexLookup.find(p);
+							if (vit != vertexLookup.end()) {
+								idc[k] = vit->second;
+							} else {
+								idc[k] = vertices.size();
+								vertices.push_back(p);
+								vertexLookup[p] = idc[k];
+							}
+						}
+						triangles.push_back(idc);
+					}
+				}
+			}
+		}
+
+
+	}
+
+};

src/TextureMapper.h

@@ -0,0 +1,201 @@
+#ifndef SRC_TEXTUREMAPPER_H_
+#define SRC_TEXTUREMAPPER_H_
+
+#include <vector>
+#include <glm/vec3.hpp>
+
+class TextureMapper {
+public:
+	typedef std::vector<glm::uint8[3]> texture3_t;
+	typedef std::vector<uint8_t> texture1_t;
+	texture3_t diffuse;
+	texture3_t normal;
+	texture1_t roughness;
+	size_t textureSize;
+
+	TextureMapper(size_t textureSize) :
+			textureSize(textureSize), diffuse(textureSize * textureSize), normal(
+					textureSize * textureSize), roughness(
+					textureSize * textureSize) {
+
+	}
+
+	void map() {
+		cout << "Stage4: UV map" << endl;
+		TriMesh mappedMesh;
+		float texf = 1. / texture_size;
+
+		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());
+	}
+};
+
+#endif /* SRC_TEXTUREMAPPER_H_ */

src/main.cpp

@@ -0,0 +1,385 @@
+#include "Density.h"
+#include "Polygoniser.h"
+
+#include <glad/glad.h>
+#include <GLFW/glfw3.h>
+#include <glm/mat4x4.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include <imgui.h>
+#include <imgui_impl_glfw_gl3.h>
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include <iostream>
+#include <fstream>
+
+static bool glck(const char* call) {
+	int err = glGetError();
+	if (err == 0) {
+		return true;
+	} else {
+		switch (err) {
+		case GL_INVALID_VALUE:
+			printf("GL_INVALID_VALUE: %s\n", call);
+			break;
+		case GL_INVALID_ENUM:
+			printf("GL_INVALID_ENUM: %s\n", call);
+			break;
+		case GL_INVALID_OPERATION:
+			printf("GL_INVALID_OPERATION: %s\n", call);
+			break;
+		default:
+			printf("UNKNOWN Error: %s\n", call);
+		}
+		return false;
+	}
+}
+
+#define GLCK(a) {a; assert(glck(#a));}
+
+// ! https://learnopengl.com/#!PBR/Lighting
+
+static const char* vertex_shader_text = "#version 130\n"
+		"uniform mat4 MVP;\n"
+		"in vec3 vPos;\n"
+		"in vec3 vNorm;\n"
+		"out vec3 normal;\n"
+		"void main()\n"
+		"{\n"
+		"    gl_Position = MVP * vec4(vPos, 1.0);\n"
+		"    normal = normalize(vNorm);\n"
+		"}\n";
+
+static const char* fragment_shader_text = "#version 130\n"
+		"in vec3 normal;\n"
+		"out vec4 color;\n"
+		"void main()\n"
+		"{\n"
+		"    color = normal.y * vec4(1.0, 0.95, 0.9, 1.0);\n"
+		"}\n";
+static void error_callback(int error, const char* description) {
+	fprintf(stderr, "Error: %s\n", description);
+}
+static void key_callback(GLFWwindow* window, int key, int scancode, int action,
+		int mods) {
+	if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS)
+		glfwSetWindowShouldClose(window, GLFW_TRUE);
+}
+
+void checkShader(GLuint vertexShader, const std::string &name) {
+	GLint status;
+	glGetShaderiv(vertexShader, GL_COMPILE_STATUS, &status);
+
+	if (status == GL_FALSE) {
+
+		std::cout << "Shader '" << name << "' compiled with an error(s):"
+				<< std::endl;
+		char buffer[512];
+		int l;
+		glGetShaderInfoLog(vertexShader, 512, &l, buffer);
+		buffer[l] = 0;
+		std::cout << buffer << std::endl;
+	}
+}
+
+void checkProgam(GLuint vertexShader, const std::string &name) {
+	GLint status;
+	glGetProgramiv(vertexShader, GL_LINK_STATUS, &status);
+
+	if (status == GL_FALSE) {
+
+		std::cout << "Program '" << name << "' failed with an error(s):"
+				<< std::endl;
+		char buffer[512];
+		int l;
+		glGetProgramInfoLog(vertexShader, 512, &l, buffer);
+		buffer[l] = 0;
+		std::cout << buffer << std::endl;
+	}
+}
+
+void calculateNormals(const std::vector<glm::vec3> &vtx,
+		std::vector<glm::uvec3> &tri, std::vector<glm::vec3> &normals) {
+	normals.clear();
+	normals.resize(vtx.size(), glm::vec3(0));
+	for (size_t i = 0; i < tri.size(); i++) {
+		int a = tri[i].x, b = tri[i].y, c = tri[i].z;
+		glm::vec3 faceNormal = normalize(
+				glm::cross(vtx[b] - vtx[a], vtx[c] - vtx[a]));
+		normals[a] += faceNormal;
+		normals[b] += faceNormal;
+		normals[c] += faceNormal;
+	}
+	for (size_t i = 0; i < normals.size(); i++)
+		normals[i] = glm::normalize(normals[i]);
+
+}
+
+void moveToMean(std::vector<glm::vec3> &vtx) {
+	glm::vec3 mean(0);
+	for (size_t i = 0; i < vtx.size(); i++) {
+		mean += vtx[i];
+	}
+	mean *= glm::vec3(1.f / vtx.size());
+	for (size_t i = 0; i < vtx.size(); i++) {
+		vtx[i] -= mean;
+	}
+}
+
+void smooth(std::vector<glm::vec3> &vtx, const std::vector<glm::uvec3> &tri) {
+	std::vector<glm::vec3> cogs(vtx.size(), glm::vec3(0.f));
+	std::vector<int> valence(vtx.size(), 0);
+
+	for (size_t iTri = 0; iTri < tri.size(); iTri++) {
+		const glm::uvec3 &idx = tri[iTri];
+		for (size_t iE = 0; iE < 3; iE++) {
+			valence[idx[iE]] += 2;
+			cogs[idx[iE]] += vtx[idx[(iE + 1) % 3]];
+			cogs[idx[iE]] += vtx[idx[(iE + 2) % 3]];
+		}
+	}
+	/*
+	 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);
+
+	 */
+	for (size_t i = 0; i < vtx.size(); i++) {
+//		vtx[i] = vtx[i] * glm::vec3(0.8)
+//				+ cogs[i] * glm::vec3(0.2f / valence[i]);
+		vtx[i] = cogs[i] * glm::vec3(1.f / valence[i]);
+	}
+}
+
+void saveAttrib(std::ostream &out, const char *prefix,
+		std::vector<glm::vec3> &elements) {
+	for (size_t i = 0; i < elements.size(); i++) {
+		out << prefix << elements[i].x << " " << elements[i].y << " " << elements[i].z
+				<< "\n";
+	}
+}
+
+void saveAttrib(std::ostream &out, const char *prefix,
+		std::vector<glm::vec2> &elements) {
+	for (size_t i = 0; i < elements.size(); i++) {
+		out << prefix << elements[i].x << " " << elements[i].y << "\n";
+	}
+}
+
+void saveFaces(std::ostream &out, const std::vector<glm::uvec3> &tris, size_t attribs) {
+	for (size_t i = 0; i < tris.size(); i++) {
+		out << "f";
+
+		for (size_t j = 0; j < 3; j++) {
+			int v = tris[i][j] + 1;
+			out << " " << v;
+			if (attribs > 1)
+				out << "/" << v;
+			if (attribs > 2)
+				out << "/" << v;
+			out << " ";
+		}
+
+		out << "\n";
+	}
+}
+
+void calculateEdges(const std::vector<glm::vec3> &vtx,
+		std::vector<glm::uvec3> &tri, std::vector<glm::vec4> &edges) {
+	edges.clear();
+}
+
+int main(void) {
+
+	std::srand(123);
+
+	Density density;
+	density.addRandomImpacts(1000, 0.2, 0.01, 0.005, 0.2, 4096);
+	//density.saveCrossSection("density.png", 512);
+
+	Polygonizer polygonizer(density);
+	polygonizer.polygonize(0.2, 1. / 20.);
+
+	moveToMean(polygonizer.vertices);
+	{
+		std::ofstream out("stage1.obj");
+		saveAttrib(out, "v ", polygonizer.vertices);
+		saveFaces(out, polygonizer.triangles, 1);
+	}
+	smooth(polygonizer.vertices, polygonizer.triangles);
+	{
+		std::ofstream out("stage2.obj");
+		saveAttrib(out, "v ", polygonizer.vertices);
+		saveFaces(out, polygonizer.triangles, 1);
+	}
+	std::vector<glm::vec3> normals;
+	calculateNormals(polygonizer.vertices, polygonizer.triangles, normals);
+	{
+		std::ofstream out("stage3.obj");
+		saveAttrib(out, "v ", polygonizer.vertices);
+		saveAttrib(out, "vn ", normals);
+		saveFaces(out, polygonizer.triangles, 2);
+	}
+
+	std::vector<glm::uvec4> edges;
+	std::cout << "Vertices: " << polygonizer.vertices.size() << std::endl;
+	std::cout << "Triangles: " << polygonizer.triangles.size() << std::endl;
+
+	GLFWwindow* window;
+	GLuint vertex_array, vertex_buffer, index_buffer, normal_buffer,
+			vertex_shader, fragment_shader, program;
+	GLint mvp_location, vpos_location, vnorm_location;
+
+	glfwSetErrorCallback(error_callback);
+
+	if (!glfwInit())
+		exit(EXIT_FAILURE);
+
+//	glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
+	//glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_COMPAT_PROFILE);
+	glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
+	glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 1);
+
+	window = glfwCreateWindow(640, 480, "Simple example", NULL, NULL);
+	if (!window) {
+		glfwTerminate();
+		exit(EXIT_FAILURE);
+	}
+
+	glfwSetKeyCallback(window, key_callback);
+
+	glfwMakeContextCurrent(window);
+	gladLoadGLLoader((GLADloadproc) glfwGetProcAddress);
+	glfwSwapInterval(1);
+
+	std::cout << "OpenGL version: " << glGetString(GL_VERSION) << std::endl;
+	std::cout << "GLSL version: " << glGetString(GL_SHADING_LANGUAGE_VERSION)
+			<< std::endl;
+	std::cout << "Vendor: " << glGetString(GL_VENDOR) << std::endl;
+	std::cout << "Renderer: " << glGetString(GL_RENDERER) << std::endl;
+
+	// NOTE: OpenGL error checks have been omitted for brevity
+
+//	    glGenBuffers(1, &vertex_buffer);
+//	    glBindBuffer(GL_ARRAY_BUFFER, vertex_buffer);
+//	    glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
+
+	GLCK(glGenVertexArrays(1, &vertex_array))
+	GLCK(glBindVertexArray(vertex_array))
+
+	GLCK(glGenBuffers(1, &vertex_buffer))
+	GLCK(glBindBuffer(GL_ARRAY_BUFFER, vertex_buffer))
+	GLCK(
+			glBufferData(GL_ARRAY_BUFFER, polygonizer.vertices.size() * sizeof(glm::vec3), polygonizer.vertices.data(), GL_STATIC_DRAW))
+
+	GLCK(glGenBuffers(1, &index_buffer))
+	GLCK(glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, index_buffer))
+	GLCK(
+			glBufferData(GL_ELEMENT_ARRAY_BUFFER, polygonizer.triangles.size() * sizeof(glm::uvec3), polygonizer.triangles.data(), GL_STATIC_DRAW))
+
+	GLCK(glGenBuffers(1, &normal_buffer))
+	GLCK(glBindBuffer(GL_ARRAY_BUFFER, normal_buffer))
+	GLCK(
+			glBufferData(GL_ARRAY_BUFFER, normals.size() * sizeof(glm::vec3), normals.data(), GL_STATIC_DRAW))
+
+	GLCK(vertex_shader = glCreateShader(GL_VERTEX_SHADER))
+	GLCK(glShaderSource(vertex_shader, 1, &vertex_shader_text, NULL))
+	GLCK(glCompileShader(vertex_shader))
+	checkShader(vertex_shader, "vertex");
+
+	GLCK(fragment_shader = glCreateShader(GL_FRAGMENT_SHADER))
+	GLCK(glShaderSource(fragment_shader, 1, &fragment_shader_text, NULL))
+	GLCK(glCompileShader(fragment_shader))
+	checkShader(fragment_shader, "fragment");
+
+	GLCK(program = glCreateProgram())
+	GLCK(glAttachShader(program, vertex_shader))
+	GLCK(glAttachShader(program, fragment_shader))
+	GLCK(glLinkProgram(program))
+	checkProgam(program, "program");
+
+	GLCK(mvp_location = glGetUniformLocation(program, "MVP"))
+	std::cout << mvp_location << std::endl;
+
+	GLCK(vpos_location = glGetAttribLocation(program, "vPos"))
+	std::cout << vpos_location << std::endl;
+
+	GLCK(vnorm_location = glGetAttribLocation(program, "vNorm"))
+	std::cout << vnorm_location << std::endl;
+
+	GLCK(glBindBuffer(GL_ARRAY_BUFFER, vertex_buffer))
+	glEnableVertexAttribArray(vpos_location);
+	GLCK(
+			glVertexAttribPointer(vpos_location, 3, GL_FLOAT, GL_FALSE, 0, (void*) 0))
+
+	if (vnorm_location >= 0) {
+		GLCK(glBindBuffer(GL_ARRAY_BUFFER, normal_buffer))
+		GLCK(glEnableVertexAttribArray(vnorm_location))
+		GLCK(
+				glVertexAttribPointer(vnorm_location, 3, GL_FLOAT, GL_FALSE, 0, (void*) 0))
+	}
+
+	ImGui_ImplGlfwGL3_Init(window, true);
+
+	while (!glfwWindowShouldClose(window)) {
+		glfwPollEvents();
+		ImGui_ImplGlfwGL3_NewFrame();
+		{
+			static float f = 0.0f;
+			static float resolution = 0.1f;
+			ImGui::Text("Hello, world!");
+			if (ImGui::Button("Wireframe"))
+				glPolygonMode( GL_FRONT_AND_BACK, GL_LINE);
+			if (ImGui::Button("Fill"))
+				glPolygonMode( GL_FRONT_AND_BACK, GL_FILL);
+			ImGui::InputFloat("Resolution", &resolution);
+			ImGui::Text("Application average %.3f ms/frame (%.1f FPS)",
+					1000.0f / ImGui::GetIO().Framerate,
+					ImGui::GetIO().Framerate);
+		}
+		int width, height;
+		glm::mat4 m(1.0f), p, mvp(2);
+
+		glfwGetFramebufferSize(window, &width, &height);
+		float ratio = width / (float) height;
+
+		glViewport(0, 0, width, height);
+		glClear(GL_COLOR_BUFFER_BIT);
+
+		m = glm::rotate(m, (float) glfwGetTime(), glm::vec3(0.f, 0.f, 1.f));
+		//m = glm::translate(m, glm::vec3(-1.f));
+		//glm::translate(mvp, glm::vec3((float) glfwGetTime()));
+		p = glm::ortho(-1.f, 1.f, -1.f, 1.f);
+		mvp = p * m;
+
+		GLCK(glUseProgram(program))
+		GLCK(
+				glUniformMatrix4fv(mvp_location, 1, GL_FALSE, (const GLfloat*) &mvp))
+
+		// Draw the triangles !
+		GLCK(glBindVertexArray(vertex_array))
+		glDrawElements( GL_TRIANGLES, polygonizer.triangles.size() * 3,
+				GL_UNSIGNED_INT, (void*) 0);
+		//GLCK(glDrawArrays( GL_TRIANGLES, 0, polygonizer.triangles.size() * 3))
+		ImGui::Render();
+		glfwSwapBuffers(window);
+
+	}
+
+	glfwDestroyWindow(window);
+
+	ImGui_ImplGlfwGL3_Shutdown();
+
+	glfwTerminate();
+}