Source

Mesh Simplification / MeshSimplification / decimation.cpp

#include "decimation.h"

//============================================================

void save_last_boundary(vector<HE_vert*>& boundary) {
	int sz;
	sz = boundary.size();
	last_boundary.clear();
	for (int i = 0 ; i < sz; i++) last_boundary.push_back(*boundary[i]);
}

void save_last_decimation(vector<HE_vert*>& boundary) {
	int sz;
	sz = boundary.size();
	last_decimated.clear();
	for (int i = 0 ; i < sz; i++) last_decimated.push_back(*boundary[i]);
}

//============================================================

int decimate(HE_vert* vert, map<int, HE_vert*>& del_verts, bool save_boundary) {
	vector<HE_vert*> boundary;
	vector<HE_face*> new_faces;
	vector<pair<pair<int, int>, int> > new_faces_points;
	int ret;
	int error_count = 0;

	ret = construct_boundary(vert, boundary);
	if (save_boundary) {
		save_last_decimation(boundary);
	}
	ret = delete_vertex(vert, del_verts);
	ret = retriangulate_boundary(vert, boundary, &new_faces, &new_faces_points);
	
	vector<HE_face*>::iterator nfiter = new_faces.begin();
	vector<pair<pair<int, int>, int> >::iterator nfpiter = new_faces_points.begin();
	for (; nfiter != new_faces.end() && nfpiter != new_faces_points.end(); nfiter++, nfpiter++) {
		if (!test_face(*nfiter, nfpiter->first.first, nfpiter->first.second, nfpiter->second)) {
			error_count--;
		}
	}

	if (error_count < 0) {
		return error_count;
	}

	return ret;
}

//============================================================

int delete_face(HE_face* f) {
	HE_edge* e = f->edge;
	HE_edge* en = e->next;

	for (int i = 0; i < 3; i++) {
		untie_edge(e, UVRT | UPAR); // untie from verts and pairs
		dele(e);
		e = en;
		en = en->next;
	}
	delf(f);

	return 0;
}

//============================================================

int delete_vertex(HE_vert* vert, map<int, HE_vert*>& del_verts) {
	HE_edge *e = vert->edge, *e_tmp;
	static map<HE_face*, int> del_faces;
	map<HE_face*, int>::iterator it;

	// delete faces around the vertex
	del_faces[e->face] = 1;
	for (e_tmp = e->next->pair; e_tmp!=e; e_tmp = e_tmp->next->pair) {
		if (e_tmp->vert != vert) {
			throw ITERATION_ERROR;
		}

		del_faces[e_tmp->face] = 1;
	}

	it = del_faces.begin();
	for (;it != del_faces.end(); it++) {
		it->second = !delete_face(it->first);
	}
	del_faces.clear();

	// delete vertex
	del_verts[vert->index] = vert;

	return 1;
}

//============================================================

int construct_boundary(HE_vert* vert, vector<HE_vert*>& boundary) {
	HE_edge *e = vert->edge, *e_tmp;

	if (e->vert != vert) {
		return ITERATION_ERROR;
	}

	// find shortest edge

	// construct boundary
	e_tmp = e->next->pair;

	boundary.push_back(e->vert2);
	for (; e_tmp != e; e_tmp = e_tmp->next->pair) {
		boundary.push_back(e_tmp->vert2);
	}

	return boundary.size();
}

//============================================================

int retriangulate_boundary(HE_vert* vert, vector<HE_vert*>& boundary, 
	vector<HE_face*>* new_faces, vector<pair<pair<int, int>, int> >* new_faces_points) {

	pair<int, int> e1, e2, e3;
	int sz = boundary.size();
	int last = sz - 1, s_last = sz - 2;
	HE_face* f;

	if (sz < 3) {
		return 0;
	}

	// retriangulate one triangle
	f = create_face(boundary[0]->index, boundary[last]->index, boundary[s_last]->index);
	if (new_faces != NULL) {
		new_faces->push_back(f);
	}
	if (new_faces_points != NULL) {
		pair<int, int> p;
		pair<pair<int, int>, int> pp;
		p.first = boundary[0]->index;
		p.second = boundary[last]->index;
		pp.first = p;
		pp.second = boundary[s_last]->index;
		new_faces_points->push_back(pp);
	}
	boundary.pop_back();

	if (sz == 3) {
		return 1;
	} else {
		return retriangulate_boundary(vert, boundary) + 1;
	}
}

//============================================================

double local_curvature(HE_vert* vert, int type) {
	double val = 0.0;

	switch (type) {
	case NORMAL_DOT_PRODUCT:
		val = normal_dot_product(vert);
		break;
	case AVERAGE_POINT:
		val = dist_to_avg_pt(vert);
		break;
	default:
		val = normal_dot_product(vert);
		break;
	}

	return val;
}

//============================================================

double normal_dot_product(HE_vert* vert) {
	HE_edge *e = vert->edge;
	Vector3 pv, p1, p2, n1, n2;
	double i = 0.0, small_len_2 = 0.0;

	pv = getVector(vert);
	p1 = getVector(e->pair->vert);
	p2 = getVector(e->next->vert);
	e = e->next->pair;
	n1 = (p2^p1).unit();

	do {
		if (e->pair == NULL ||
			e->pair->next == NULL) {
				throw ITERATION_ERROR;
		}
		p1 = p2;
		p2 = getVector(e->next->vert);
		e = e->next->pair;
		n2 = (p2^p1).unit();
		small_len_2 = n1*n2;
		n1 = n2;
	} while (e != vert->edge);

	return small_len_2;
}

//============================================================
double dist_to_avg_pt(HE_vert* vert) {
	HE_edge *e = vert->edge;
	Vector3 v, vp, p, avgp(0.0, 0.0, 0.0);
	double i = 0.0, avg_len = 0.0, curviture;

	vp = getVector(vert);

	do {
		p = getVector(e->pair->vert);
		avgp += p;
		avg_len += (p - vp).length2();
		i+=1.0;
		if (e->pair == NULL ||
			e->pair->next == NULL) {
				throw ITERATION_ERROR;
		}
		e = e->pair->next->next;
	} while (e != vert->edge);
	
	avg_len /= i;
	avgp *= (1.0/i);
	v = vp - avgp;

	curviture = v.length2();

	return (curviture*curviture)/avg_len;
}

//============================================================
map<int, HE_vert*> del_verts;
void clean_decimation_structures() {
	map<int, HE_vert*>::iterator j = del_verts.begin();
	for (; j != del_verts.end(); j++) {
		delv(j->second);
	}
}

int decimate(map<int, HE_vert*>::iterator& i, double threshold, int test_type, bool save_boundary) {
	if (decimation_test(i, threshold, test_type)) {
		decimate(i->second, del_verts, save_boundary);
		return 1;
	}

	return 0;
}

int decimation_test(map<int, HE_vert*>::iterator& i, double threshold, int test_type) {
	double test = 0, tc2 = threshold;
	tc2*=tc2;

	test = local_curvature(i->second, test_type);
	
	return test < tc2;
}

void decimate(double threshold, int test_type, bool save_boundary) {
	int total = 0, successful = 0, failed = 0, invalid1 = 0, invalid2 = 0;
	double res = 0;
	
	map<int, HE_vert*>::iterator i = verts.begin();

	cout << "Applying vertex decimation algorithm... " << endl;
	for (; i != verts.end(); i++) {
		if (del_verts.count(i->first)) continue;
		try {
			total += decimate(i, threshold, test_type, save_boundary);
		} catch (int i) {
			if (i == ITERATION_ERROR) {
				cout << "Iteration error occured in local_curvature() function." << endl;
				break;
			}
		}
	}

	clean_decimation_structures();

	cout << "Vertices decimated " << total << " of " << nv << endl;
}