meshsimplify.cpp

struct Quadric {
	float q[10];
	Quadric() = default;
	Quadric(float q0, float q1, float q2, float q3, float q4,
		float q5, float q6, float q7, float q8, float q9)
	{
		q[0] = q0;
		q[1] = q1;
		q[2] = q2;
		q[3] = q3;
		q[4] = q4;
		q[5] = q5;
		q[6] = q6;
		q[7] = q7;
		q[8] = q8;
		q[9] = q9;
	}

	Quadric &operator+=(const Quadric &rhs)
	{
		for (int i = 0; i < 10; i++)
			q[i] += rhs[i];
		return *this;
	}

	float &operator[](int i) { return q[i]; }
	float operator[](int i) const { return q[i]; }
};

static inline bool operator==(const Quadric &lhs, const Quadric &rhs)
{
	for (int i = 0; i < 10; i++) {
		if (fabs(lhs[0] - rhs[0]) > MATH_EPSILON)
			return false;
	}
	return true;
}

static inline bool operator!=(const Quadric &lhs, const Quadric &rhs)
{
	return !operator==(lhs, rhs);
}

static inline Quadric operator+(const Quadric &lhs, const Quadric &rhs)
{
	Quadric out;
	for (int i = 0; i < 10; i++)
		out[i] = lhs[i] + rhs[i];
	return out;
}

static inline Quadric Quadric_Zero()
{
	Quadric out;
	for (int i = 0; i < 10; i++)
		out[i] = 0;
	return out;
}

static inline Mat4 ToMat4(const Quadric &q)
{
	return Mat4(
		q[0], q[1], q[2], q[3],
		q[1], q[4], q[5], q[6],
		q[2], q[5], q[7], q[8],
		q[3], q[6], q[8], q[9]
	);
}

static inline Quadric PlaneToQuadric(const Vec4 &p)
{
	float a2 = p.x * p.x;
	float b2 = p.y * p.y;
	float c2 = p.z * p.z;
	float d2 = p.w * p.w;
	float ab = p.x * p.y;
	float ac = p.x * p.z;
	float ad = p.x * p.w;
	float bc = p.y * p.z;
	float bd = p.y * p.w;
	float cd = p.z * p.w;
	return Quadric(a2, ab, ac, ad, b2, bc, bd, c2, cd, d2);
}

struct HalfEdge {
	int vertex;
	int face;
	int next;
};

enum VertexFlags {
	VF_VIRTUAL = 1 << 0,
	VF_PINNED = 1 << 1,
	VF_BOUNDARY = 1 << 2,
};

struct V3M1_HE {
	Vec3 position;
	uint16_t material;
	uint16_t flags;
	int he;
	Quadric quadric;
};

struct Face {
	int he;
	Vec4 plane;
};

struct Edge {
	int start;
	int end;
};

static inline bool operator==(const Edge &lhs, const Edge &rhs)
{
	return lhs.start == rhs.start && lhs.end == rhs.end;
}
static inline bool operator!=(const Edge &lhs, const Edge &rhs)
{
	return lhs.start != rhs.start || lhs.end != rhs.end;
}

int Hash(const Edge &e)
{
	return Hash(Slice<const char>((const char*)&e, sizeof(e)));
}

struct RankedEdge {
	Vec3 new_position;
	float rank;
	int he;
};

struct Triangle {
	Vec3 a;
	Vec3 b;
	Vec3 c;
};

static inline uint32_t FastRand(uint32_t *state)
{
	uint32_t x = *state;
	x += x;
	if (x & 0x80000000L)
		x ^= 0x88888EEFUL;
	*state = x;
	return x;
}

constexpr float CANT_COLLAPSE_RANK = 100.0f;

struct HalfEdgeMesh {
	Vector<HalfEdge> m_edges;
	Vector<V3M1_HE> m_vertices;
	Vector<Face> m_faces;
	HashMap<Edge, int> m_edge_map;
	Vector<Triangle> m_triangles;

	void Dump()
	{
		printf("Edges size: %d bytes (%d)\n", m_edges.ByteLength(), m_edges.Length());
		printf("Vertices size: %d bytes (%d)\n", m_vertices.ByteLength(), m_vertices.Length());
		printf("Faces size: %d bytes (%d)\n", m_faces.ByteLength(), m_faces.Length());
	}

	HalfEdgeMesh &Init()
	{
		m_edges.Init();
		m_vertices.Init();
		m_faces.Init();
		m_edge_map.Init();
		m_triangles.Init();
		return *this;
	}

	void Free()
	{
		m_edges.Free();
		m_vertices.Free();
		m_faces.Free();
		m_edge_map.Free();
		m_triangles.Free();
	}

	/*
	HalfEdge &HE(int n) { return m_edges[n]; }
	Face &F(int n) { return m_faces[n]; }
	V3M1_HE &V(int n) { return m_vertices[n]; }
	*/
	HalfEdge &HE(int n) { return m_edges.Data()[n]; }
	Face &F(int n) { return m_faces.Data()[n]; }
	V3M1_HE &V(int n) { return m_vertices.Data()[n]; }

	void DebugDrawFace(int n, const Vec3 &off, const Vec3 &color = Vec3_X())
	{
		if (n == -1) {
			Warn("drawing non-existent face");
			return;
		}
		const int h0 = m_faces[n].he;
		const int h1 = HE(h0).next;
		const int h2 = HE(h1).next;
		const Vec3 v0 = V(HE(h0).vertex).position;
		const Vec3 v1 = V(HE(h1).vertex).position;
		const Vec3 v2 = V(HE(h2).vertex).position;
		printf("Drawing face at: (%f %f %f), (%f %f %f), (%f %f %f)\n", VEC3(v0), VEC3(v1), VEC3(v2));
		debug_draw.Line(v0+off, v1+off, color);
		debug_draw.Line(v1+off, v2+off, color);
		debug_draw.Line(v2+off, v0+off, color);
	}

	void DebugDrawEdge(int n, const Vec3 &off, const Vec3 &color = Vec3_X())
	{
		printf("Drawing edge %d\n", n);
		const int h0 = n;
		const int h1 = h0^1;
		const Vec3 v0 = V(HE(h0).vertex).position;
		const Vec3 v1 = V(HE(h1).vertex).position;
		debug_draw.Line(v0+off, v1+off, color);
	}

	void UpdateVertexQuadric(int n)
	{
		V3M1_HE &v = V(n);
		v.quadric = EdgeVertexQuadric(v.he);
	}

	void UpdateNeighbourVertexQuadrics(int n)
	{
		const int h0 = V(n).he;
		int he = h0;
		do {
			UpdateVertexQuadric(HE(he^1).vertex);
			he = HE(he).next^1;
		} while (he != h0);
	}

	void UpdateFacePlane(int n)
	{
		Face &f = F(n);
		const int a = f.he;
		const int b = HE(a).next;
		const int c = HE(b).next;
		const Vec3 va = V(HE(a).vertex).position;
		const Vec3 vb = V(HE(b).vertex).position;
		const Vec3 vc = V(HE(c).vertex).position;
		const Vec3 ab = va - vb;
		const Vec3 cb = vc - vb;
		const Vec3 N = Normalize(Cross(cb, ab));
		const float d = -Dot(N, va);
		f.plane = Vec4(N.x, N.y, N.z, d);
	}

	void UpdateNeighbourFacePlanes(int n)
	{
		const int h0 = V(n).he;
		int he = h0;
		do {
			UpdateFacePlane(HE(he).face);
			he = HE(he).next^1;
		} while (he != h0);
	}

	void CollapseRandomEdges(float threshold)
	{
		constexpr int N_EDGES = 5;
		uint32_t rnd = 0x49f6428aUL;
		int n_failed = 0;

		if (threshold >= CANT_COLLAPSE_RANK)
			threshold = CANT_COLLAPSE_RANK - 1.0f;

		for (int i = 0; i < m_faces.Length(); i++)
			UpdateFacePlane(i);
		for (int i = 0; i < m_vertices.Length(); i++) {
			if (IsBoundaryEdgeVertex(V(i).he))
				V(i).flags |= VF_BOUNDARY;
			UpdateVertexQuadric(i);
		}

		while (n_failed < 100) {
			RankedEdge the_edge;
			the_edge.rank = CANT_COLLAPSE_RANK;
			const uint32_t rnd_div = m_edges.Length()/2;
			for (int i = 0; i < N_EDGES; i++) {
				int j = (FastRand(&rnd) % rnd_div)*2;
				const RankedEdge e = RankEdge(j);
				if (e.rank < the_edge.rank)
					the_edge = e;
			}

			if (the_edge.rank < threshold) {
				CollapseEdge(the_edge.he, the_edge.new_position);
				n_failed = 0;
			} else {
				n_failed++;
			}
		}
	}

	RankedEdge RankEdge(int n)
	{
		RankedEdge re;
		if (!CanCollapseEdge(n)) {
			re.rank = CANT_COLLAPSE_RANK;
			return re;
		}

		const int h0 = n;
		const int h1 = h0^1;
		const Mat4 Q = ToMat4(V(HE(h0).vertex).quadric + V(HE(h1).vertex).quadric);
		Mat4 Qtmp = Q;
		Qtmp[3]  = 0;
		Qtmp[7]  = 0;
		Qtmp[11] = 0;
		Qtmp[15] = 1;

		bool ok = false;
		Mat4 Qi = Inverse(Qtmp, &ok);
		if (ok) {
			Vec4 p = Qi * Vec4(0, 0, 0, 1);
			re.new_position = ToVec3(p);
			if (IsNaN(re.new_position)) {
				re.new_position = EdgeMiddlePoint(n);
				p = ToVec4(re.new_position);
			}
			re.rank = fabs(Dot(p * Q, p));
			re.he = n;
		} else {
			const Vec4 mp = ToVec4(EdgeMiddlePoint(n));
			re.new_position = ToVec3(mp);
			re.rank = fabs(Dot(mp * Q, mp));
			re.he = n;
		}

		if (!CanCollapseEdgeAt(n, re.new_position))
			re.rank += CANT_COLLAPSE_RANK;
		return re;
	}

	Vec3 EdgeMiddlePoint(int n)
	{
		const int h0 = n;
		const int h1 = h0^1;
		const Vec3 v0 = V(HE(h0).vertex).position;
		const Vec3 v1 = V(HE(h1).vertex).position;
		return (v0 + v1) / Vec3(2);
	}

	Quadric FaceQuadric(int n)
	{
		return PlaneToQuadric(F(n).plane);
	}

	Quadric EdgeVertexQuadric(int n)
	{
		Quadric q = Quadric_Zero();
		const int h0 = n;
		int he = h0;
		bool open = false;
		do {
			if (HE(he).face == -1) {
				open = true;
				break;
			} else {
				q += FaceQuadric(HE(he).face);
			}
			he = HE(he).next^1;
		} while (he != h0);
		if (open) {
			he = h0^1;
			while (HE(he).face != -1) {
				q += FaceQuadric(HE(he).face);
				he = HE(HE(he).next).next^1;
			}
		}
		return q;
	}

	void PostCollapse(int v1)
	{
		UpdateVertexQuadric(v1);
		UpdateNeighbourVertexQuadrics(v1);
		UpdateNeighbourFacePlanes(v1);
	}

	void CollapseEdge3Triangles(int n, const Vec3 &pos)
	{
		const int h0 = n;
		const int h1 = n^1;
		const int v1 = HE(h1).vertex;
		const int a = HE(h0).next;
		const int b = HE(h1).next;
		const int c = HE(b).next;
		const int d = HE(a).next;
		const int e = c^1;
		const int f = a^1;
		const int g = HE(e).next;

		// move h1's vertex to new position
		V(v1).position = pos;

		// redirect all h0's vertex edges to h1's vertex
		HE(c).vertex = v1;
		HE(f).vertex = v1;

		// fix remaining face
		F(HE(g).face).he = g;

		// reconnect edges
		HE(g).next = d;
		HE(d).next = b;
		HE(b).next = g;

		// reconnect face
		HE(b).face = HE(g).face;
		HE(d).face = HE(g).face;

		// fix vertices
		V(HE(h1).vertex).he = d;
		V(HE(e).vertex).he = b;
		V(HE(a).vertex).he = g;

		// erase main and side edges
		int rm[] = {h0, c, a};
		EraseEdges(rm);
		PostCollapse(v1);
	}

	void RedirectEdgesTo(int edge, int newvertex)
	{
		int he = edge;
		do {
			HE(he).vertex = newvertex;
			he = HE(he).next^1;
		} while (he != edge);
	}

	bool IsBoundaryEdgeVertex(int n)
	{
		bool open = false;
		int he = n;
		do {
			if (HE(he).face == -1) {
				open = true;
				break;
			}
			he = HE(he).next^1;
		} while (he != n);
		return open;
	}

	void EraseEdge(int n)
	{
		int cur = n & ~1;
		int last = (m_edges.Length() - 1) & ~1;
		if (last != cur)
			MoveEdge(last, cur);
		m_edges.Resize(m_edges.Length() - 2);
	}

	void EraseEdges(int *rm)
	{
		if (rm[0] < rm[1]) std::swap(rm[0], rm[1]);
		if (rm[1] < rm[2]) std::swap(rm[1], rm[2]);
		if (rm[0] < rm[1]) std::swap(rm[0], rm[1]);

		for (int i = 0; i < 3; i++)
			EraseEdge(rm[i]);
	}

	void MoveEdge(int from, int to)
	{
		const int from2 = from^1;
		const int to2 = to^1;
		const int f0 = HE(from).face;
		const int f1 = HE(from2).face;
		if (f0 != -1) {
			if (F(f0).he == from) {
				F(f0).he = to;
			}
			int he = HE(HE(from).next).next;
			HE(he).next = to;
		}
		if (f1 != -1) {
			if (F(f1).he == from2) {
				F(f1).he = to2;
			}
			int he = HE(HE(from2).next).next;
			HE(he).next = to2;
		}
		const int v0 = HE(from).vertex;
		const int v1 = HE(from2).vertex;
		if (V(v0).he == from)
			V(v0).he = to;
		if (V(v1).he == from2)
			V(v1).he = to2;

		m_edges[to] = m_edges[from];
		m_edges[to^1] = m_edges[from^1];
	}

	void MoveFace(int from, int to)
	{
		int h0 = F(from).he;
		int h1 = HE(h0).next;
		int h2 = HE(h1).next;
		HE(h0).face = to;
		HE(h1).face = to;
		HE(h2).face = to;
		m_faces[to] = m_faces[from];
	}

	void EraseFace(int n)
	{
		const int last = m_faces.Length()-1;
		if (last != n)
			MoveFace(last, n);
		m_faces.Resize(m_faces.Length()-1);
	}

	void EraseFaces(int *rm)
	{
		if (rm[0] < rm[1]) std::swap(rm[0], rm[1]);
		EraseFace(rm[0]);
		EraseFace(rm[1]);
	}

	void ReassignVertex(int from, int to)
	{
		const int h0 = V(from).he;
		int he = h0;
		bool open = false;
		do {
			HE(he).vertex = to;
			if (HE(he).face == -1) {
				open = true;
				break;
			}
			he = HE(he).next^1;
		} while (he != h0);

		if (open) {
			he = h0^1;
			while (HE(he).face != -1) {
				he = HE(HE(he).next).next;
				HE(he).vertex = to;
				he = he^1;
			}
		}
	}

	void MoveVertex(int from, int to)
	{
		ReassignVertex(from, to);
		m_vertices[to] = m_vertices[from];
	}

	void SwapVertices(int a, int b)
	{
		ReassignVertex(a, b);
		ReassignVertex(b, a);
		std::swap(m_vertices[a], m_vertices[b]);
	}

	// returns the length of vertices without virtual ones, this is also the
	// index of the first virtual vertex
	int MoveVirtualVerticesToEnd()
	{
		int last_non_virtual = -1;
		for (int i = m_vertices.Length()-1; i >= 0; i--) {
			if ((V(i).flags & VF_VIRTUAL) == 0) {
				last_non_virtual = i;
				break;
			}
		}
		if (last_non_virtual == -1)
			return 0;

		int first_virtual = -1;
		for (int i = 0; i < m_vertices.Length(); i++) {
			if (V(i).flags & VF_VIRTUAL) {
				first_virtual = i;
				break;
			}
		}
		if (first_virtual == -1)
			return m_vertices.Length();

		while (first_virtual < last_non_virtual) {
			SwapVertices(first_virtual, last_non_virtual);
			while (first_virtual < last_non_virtual && (V(first_virtual).flags & VF_VIRTUAL) == 0)
				first_virtual++;
			while (first_virtual < last_non_virtual && V(last_non_virtual).flags & VF_VIRTUAL)
				last_non_virtual--;
		}

		return first_virtual;
	}

	void EraseVertex(int n)
	{
		const int last = m_vertices.Length()-1;
		if (last != n)
			MoveVertex(last, n);
		m_vertices.Resize(m_vertices.Length()-1);
	}

	void CollectTriangle(int n)
	{
		const int h0 = n;
		const int h1 = HE(h0).next;
		const int h2 = HE(h1).next;
		const Vec3 a = V(HE(h0).vertex).position;
		const Vec3 b = V(HE(h1).vertex).position;
		const Vec3 c = V(HE(h2).vertex).position;
		m_triangles.Append({a, b, c});
	}

	// collect all triangles around the n's vertex, except two triangles
	// next to the 'n' edge, 'a' vertex is always the vertex of 'n'
	void CollectVertexTriangles(int n)
	{
		int he = HE(n).next^1;
		int end = HE(HE(n^1).next).next;
		do {
			CollectTriangle(he);
			he = HE(he).next^1;
		} while (he != end);
	}

	bool Validate()
	{
		for (int i = 0; i < m_edges.Length(); i++) {
			if (HE(i).face >= m_faces.Length()) {
				Warn("%d edge's face is broken (%d)", i, HE(i).face);
				return false;
			}
		}
		return true;
	}

	bool CanCollapseEdgeAt(int n, const Vec3 &pos)
	{
		// check if normals are not flipping
		const int h0 = n;
		const int h1 = n^1;
		m_triangles.Clear();
		CollectVertexTriangles(h0);
		CollectVertexTriangles(h1);
		for (const auto &tri : m_triangles) {
			const Vec3 alt_ab = pos - tri.b;
			const Vec3 ab = tri.a - tri.b;
			const Vec3 cb = tri.c - tri.b;
			const Vec3 n0 = Cross(cb, ab);
			const Vec3 n1 = Cross(cb, alt_ab);
			if (Dot(n0, n1) < 0) {
				return false;
			}
		}

		return true;
	}

	bool CanCollapseEdge(int n)
	{
		const int h0 = n;
		const int h1 = n^1;
		if (V(HE(h0).vertex).flags & (VF_PINNED | VF_BOUNDARY))
			return false;
		if (V(HE(h1).vertex).flags & (VF_PINNED | VF_BOUNDARY))
			return false;

		/*
		const bool open0 = IsBoundaryEdgeVertex(h0);
		const bool open1 = IsBoundaryEdgeVertex(h1);

		// one of the vertices is boundary
		if (open0 || open1)
			return false;
		*/

		// avoid closing triangle holes, leads to non-manifold geometry
		{
			const int v0 = HE(HE(h0).next).vertex;
			const int v1 = HE(HE(h1).next).vertex;
			int he = h0;
			do {
				int he2 = h1;
				do {
					int v = HE(he^1).vertex;
					if (v == HE(he2^1).vertex && v != v0 && v != v1)
						return false;
					he2 = HE(he2).next^1;
				} while (he2 != h1);
				he = HE(he).next^1;
			} while (he != h0);
		}
		/*
		// XXX(nsf): seems like it's not enough, but I haven't seen a case
		// where it failed
		{
			const int a = HE(h1).next;
			const int f = HE(a).next;
			const int L = HE(a^1).next;
			const int g = HE(f^1).next;
			const int h = HE(g).next;
			if (HE(L).vertex == HE(h^1).vertex)
				return false;

			const int b = HE(h0).next;
			const int e = HE(b).next;
			const int c = HE(b^1).next;
			const int i = HE(e^1).next;
			const int k = HE(i).next;
			if (HE(c).vertex == HE(k^1).vertex)
				return false;
		}
		*/

		// these checks also handle tetrahedron cases, they are not handled
		// by triangle hole test above
		{
			const int b = HE(h0).next;
			const int e = HE(b).next;
			const int d = HE(HE(b^1).next).next;
			const int i = HE(e^1).next;
			if (i == (d^1)) {
				return false;
			}
		}
		{
			const int a = HE(h1).next;
			const int f = HE(a).next;
			const int j = HE(HE(a^1).next).next;
			const int g = HE(f^1).next;
			if (g == (j^1)) {
				return false;
			}
		}

		return true;
	}

	// Assumes CanCollapseEdge is true
	bool CollapseEdge(int n, const Vec3 &pos)
	{
		const int h0 = n;
		const int h1 = h0^1;

		EraseVertex(HE(h0).vertex);

		// erase faces
		int rm_faces[] = {HE(h0).face, HE(h1).face};
		EraseFaces(rm_faces);

		const int v1 = HE(h1).vertex;

		// this is an unique case, where tri loop around h0's vertex is closed
		// and has only 3 triangles
		{
			const int left = HE(HE(h1).next).next^1;
			const int right = HE(h0).next^1;
			if (HE(left).face != -1 && HE(left).face == HE(right).face) {
				CollapseEdge3Triangles(n, pos);
				return true;
			}
		}

		// move vertex to new position
		V(v1).position = pos;

		// redirect h0's vertex edges to h1's vertex
		RedirectEdgesTo(h0, v1);

		const int b0 = HE(h0).next;
		const int b1 = b0^1;
		const int e = HE(b0).next;

		// full case, two triangles of the tri loop on the right side
		const int c = HE(b1).next;
		const int d = HE(c).next;

		// redirect edges
		HE(d).next = e;
		HE(e).next = c;

		// redirect and fix face
		F(HE(b1).face).he = e;
		HE(e).face = HE(b1).face;

		const int a = HE(h1).next;
		const int f0 = HE(a).next;
		const int f1 = f0^1;

		// full case, two triangles of the tri loop on the left side
		const int g = HE(f1).next;
		const int h = HE(g).next;

		// redirect edges
		HE(h).next = a;
		HE(a).next = g;

		// redirect and fix face
		F(HE(f1).face).he = a;
		HE(a).face = HE(f1).face;

		// fix vertex
		V(HE(h1).vertex).he = e;
		V(HE(f1).vertex).he = a;
		V(HE(b0).vertex).he = d;

		int rm_edges[] = {b0, f0, h0};
		EraseEdges(rm_edges);
		PostCollapse(v1);
		return true;
	}

	int InsertEdgeMaybe(Edge e)
	{
		const bool opposite = e.start > e.end;
		if (opposite)
			std::swap(e.start, e.end);

		const int *i = m_edge_map.Get(e);
		if (i) {
			return opposite ? *i ^ 1 : *i;
		}

		const int first = m_edges.Length();
		HalfEdge *he0 = m_edges.Append();
		he0->face = -1;
		he0->next = -1;
		he0->vertex = e.end;
		V(e.end).he = first;

		const int second = m_edges.Length();
		HalfEdge *he1 = m_edges.Append();
		he1->face = -1;
		he1->next = -1;
		he1->vertex = e.start;
		V(e.start).he = second;

		m_edge_map.Insert(e, first);
		return opposite ? second : first;
	}

	void AddTriangle(int a, int b, int c, bool pinned = false)
	{
		const int hei0 = InsertEdgeMaybe({a, b});
		const int hei1 = InsertEdgeMaybe({b, c});
		const int hei2 = InsertEdgeMaybe({c, a});
		HalfEdge *he0 = &m_edges[hei0];
		HalfEdge *he1 = &m_edges[hei1];
		HalfEdge *he2 = &m_edges[hei2];
		if (pinned) {
			V(he0->vertex).flags |= VF_PINNED;
			V(he1->vertex).flags |= VF_PINNED;
			V(he2->vertex).flags |= VF_PINNED;
		}
		const int fi = m_faces.Length();
		Face *f = m_faces.Append();

		f->he = hei2;
		he0->face = fi;
		he1->face = fi;
		he2->face = fi;
		he0->next = hei1;
		he1->next = hei2;
		he2->next = hei0;
	}

	int AddVertex(const V3M1_HE &v)
	{
		const int index = m_vertices.Length();
		V3M1_HE *nv = m_vertices.Append();
		*nv = v;
		return index;
	}

	void Build(Slice<const V3N3M1> vertices, Slice<const uint32_t> indices)
	{
		m_edges.Clear();
		m_vertices.Clear();
		m_faces.Clear();
		m_edge_map.Clear();

		m_vertices.Resize(vertices.length);
		for (int i = 0; i < vertices.length; i++) {
			m_vertices[i].position = vertices[i].position;
			m_vertices[i].material = vertices[i].material;
		}

		for (int i = 0; i < indices.length; i += 3) {
			AddTriangle(indices[i+0], indices[i+1], indices[i+2]);
		}
		m_edge_map.Free();
		m_edge_map.Init();
	}

	void Output(Vector<V3N3M1> &vertices, Vector<uint32_t> &indices,
		Vector<Vec3> &virt_vertices, Vector<uint32_t> &virt_indices)
	{
		const int non_virt_len = MoveVirtualVerticesToEnd();
		const int virt_len = m_vertices.Length() - non_virt_len;

		vertices.Resize(non_virt_len);
		for (int i = 0; i < non_virt_len; i++) {
			vertices[i].position = m_vertices[i].position;
			vertices[i].material = m_vertices[i].material;
			vertices[i].normal = Vec3(0);
		}
		virt_vertices.Resize(virt_len+1);
		for (int i = 0; i < virt_len; i++)
			virt_vertices[1+i] = m_vertices[non_virt_len + i].position;

		indices.Reserve(m_faces.Length() * 3);
		indices.Clear();
		virt_indices.Clear();
		for (const auto &face : m_faces) {
			int he = face.he;
			const int v0 = HE(he).vertex;
			he = HE(he).next;
			const int v1 = HE(he).vertex;
			he = HE(he).next;
			const int v2 = HE(he).vertex;
			const bool v0_is_virt = v0 >= non_virt_len;
			const bool v1_is_virt = v1 >= non_virt_len;
			const bool v2_is_virt = v2 >= non_virt_len;
			if (v0_is_virt || v1_is_virt || v2_is_virt) {
				const int v0i = v0_is_virt ? -(v0 - non_virt_len) - 1 : v0;
				const int v1i = v1_is_virt ? -(v1 - non_virt_len) - 1 : v1;
				const int v2i = v2_is_virt ? -(v2 - non_virt_len) - 1 : v2;
				virt_indices.Append(v0i);
				virt_indices.Append(v1i);
				virt_indices.Append(v2i);
			} else {
				indices.Append(v0);
				indices.Append(v1);
				indices.Append(v2);
			}
		}
	}
};

/*

	half_edge_mesh.Build(vertices0, indices0);
	printf("HalfEdgeMesh size: %d bytes\n",
		half_edge_mesh.m_edges.ByteLength() +
		half_edge_mesh.m_faces.ByteLength() +
		half_edge_mesh.m_vertices.ByteLength());

	auto &m = half_edge_mesh;

	if (threshold != 0.0f)
		m.CollapseRandomEdges(threshold * threshold);

	printf("HalfEdgeMesh size (after simplification): %d bytes\n",
		half_edge_mesh.m_edges.ByteLength() +
		half_edge_mesh.m_faces.ByteLength() +
		half_edge_mesh.m_vertices.ByteLength());

	half_edge_mesh.Output(vertices, indices);
	printf("Mesh size: %d bytes (threshold: %f)\n", vertices.ByteLength() + indices.ByteLength(), threshold);

*/