decimate.cpp

#include <ultimaille/all.h>
#include <set>
#define FOR(i, n) for(int i = 0; i < static_cast<int>(n); i++)

using namespace UM;



struct EdgeCollapse {
	Triangles& tri;
	Surface::Halfedge h;
	Surface::Halfedge opp;
	EdgeCollapse(Triangles& tri, Surface::Halfedge h) :tri(tri), h(h), opp(h.opposite()) { }

	bool valid() {
		// this version does not decimate boundaries
		for (auto cir : h.from().iter_halfedges()) if (!cir.opposite().active()) return false;
		for (auto cir : h.to().iter_halfedges()) if (!cir.opposite().active()) return false;
		
		// avoid producing edges cycles of length 2
		for (auto cir_from : h.from().iter_halfedges())
			for (auto cir_to : h.to().iter_halfedges())
				if (cir_from.to() == cir_to.to())
					if (cir_from.to() != h.next().to())
						if (!opp.active() || cir_from.to() != opp.next().to())
							return false;
		return true;
	}

	void apply() {
		int to = h.to();
		int facets_to_kill[2] = {h.facet(),opp.facet()};
		std::vector<int> c_incident;
		for (auto cir_from : h.from().iter_halfedges()) c_incident.push_back(cir_from);
		for (int c : c_incident) tri.conn->change_from(Surface::Halfedge(tri, c), to);
		for(int f: facets_to_kill) tri.conn->active[f] = false;
	}
};


namespace CollapseVolumePreserving {
	mat4x4 quadratic_form(Triangles& tri, Surface::Halfedge h) {
		std::set<int> triangles;
		for (auto cir : h.from().iter_halfedges())triangles.insert(cir.facet());
		for (auto cir : h.to().iter_halfedges())triangles.insert(cir.facet());
		mat4x4 Q = {};
		for (int tid : triangles) {
			Triangle3 tr(Surface::Facet(tri, tid));
			vec3 n = tr.normal();
			vec3 G = tr.bary_verts();
			vec4 p{ n[0], n[1], n[2], -n * G };
			FOR(i, 4)FOR(j, 4) Q[i][j] += p[i] * p[j];
		}
		return Q;
	}
	std::tuple<vec3,double> pos_and_cost(Triangles& tri, Surface::Halfedge h) {
		EdgeCollapse ec(tri, h);
		if (!ec.valid()) return { vec3(),1e20 };
		mat4x4 Q = quadratic_form(tri, h);
		mat4x4 A = Q;
		
		FOR(j, 3) A[3][j] = 0;
		A[3][3] = 1;
		mat4x4 M = A.invert();
		vec4 pos = M * vec4{ 0, 0, 0, 1 };
		vec3 P(pos[0], pos[1], pos[2]);
		double nrj = pos * (Q * pos);
		return { P,nrj };
	}

	vec3 vertex_position(Triangles& tri, Surface::Halfedge h) {
		auto [P, nrj] = pos_and_cost(tri, h);
		return P;
	};


	double cost(Triangles& tri, Surface::Halfedge h) {
		auto [P, nrj] = pos_and_cost(tri, h);
		return nrj;
	};
};



namespace CollapseSmallEdges {
	vec3 vertex_position([[maybe_unused]] Triangles& tri, Surface::Halfedge h) {
		return h.to().pos();
	};

	double cost(Triangles& tri, Surface::Halfedge h) {
		EdgeCollapse ec(tri, h);
		if (!ec.valid()) return 1e20;
		return Segment3(h).length();
	};
}



// simplest example
//using namespace CollapseSmallEdges;

// using quadratic 
//using namespace CollapseVolumePreserving;


bool use_quadratic_objective = true;
vec3 vertex_position(Triangles& tri, Surface::Halfedge h) {
	if (use_quadratic_objective) return CollapseVolumePreserving::vertex_position(tri, h);
	else return CollapseSmallEdges::vertex_position(tri, h);
};

double cost(Triangles& tri, Surface::Halfedge h) {
	if (use_quadratic_objective) return CollapseVolumePreserving::cost(tri, h);
	else return CollapseSmallEdges::cost(tri, h);
};


void decimate(Triangles& tri, PointAttribute<bool>& selection, double reduction_factor = 10) {
	int nb_max_collapse = static_cast<int>((1. - 1. / reduction_factor) * tri.nverts());

	CornerAttribute<double> stored_cost(tri, 1e20);
	for (auto h : tri.iter_halfedges()) {
		if (!selection[h.from()] && !selection[h.to()]) {
			stored_cost[h] = cost(tri,h);
		} else {
			std::cout << "avoided\n";
		}
	};


	auto cmp = [&](int i, int j) {
		if (i == -1) return true;
		if (j == -1) return false;
		return stored_cost[i] < stored_cost[j];
	};

	std::set<int, decltype(cmp)> heap(cmp);
	for (auto h : tri.iter_halfedges()) if (stored_cost[h] < 1e20) heap.insert(h);

	while (!heap.empty()) {
		auto cur = heap.upper_bound(-1);
		Surface::Halfedge h(tri, *cur);

		heap.erase(cur);
		if (!h.active()) continue;
		if (nb_max_collapse-- < 0) break;

		EdgeCollapse ec(tri, h);
		h.to().pos() = vertex_position(tri,h);
		auto v = h.to();
		ec.apply();

		
		std::set<int> h_to_update;
		// find the set of halfedge with modified cost
		for (auto iter_h : v.iter_halfedges()) for (auto cir : iter_h.to().iter_halfedges()) {
			h_to_update.insert(cir);
			auto opp = cir.opposite();
			if (opp.active()) h_to_update.insert(opp);
		}

		// update the cost and the heap accordingly
		for (int it : h_to_update) {
			auto handle = heap.find(it);
			if (handle != heap.end()) heap.erase(handle);
			stored_cost[it] = cost(tri, Surface::Halfedge(tri,it));
			if (stored_cost[it] < 1e20)
				heap.insert(it);
		}
	}

	tri.compact();
}

void shift() {
    UM::Triangles m;
    SurfaceAttributes attributes = read_by_extension("../pre.geogram", m);
    PointAttribute<bool> selection("selection", attributes, m);
    m.connect();

    for (auto v : m.iter_vertices()) {
        if (selection[v]) {
            m.points[v].x += 0.25;
        }
    }

    for (auto i = 0; i < 1000; i++) {
        for (auto v: m.iter_vertices()) {
            if (!selection[v] && !v.on_boundary()) {
                vec4 bary;
                for (auto he : v.iter_halfedges()) {
                    bary += he.to().pos().xyz1();
                }
                m.points[v] = bary.xyz() / bary[3];
            }
        }
    }
    write_by_extension("after.geogram", m);
}


int main() {
	FOR(iter, 2) {
		use_quadratic_objective = (iter == 1);
		shift();
		Triangles m;
		Triangles m1;
		SurfaceAttributes attributes = read_by_extension("../pre.geogram", m1);
		PointAttribute<bool> selection("selection", attributes, m1);
		read_by_extension("after.geogram", m);
		m.connect();
		decimate(m, selection, 1.1);
        write_by_extension("after_after" + std::to_string(iter) + ".geogram", m);
	}
    return 0;
}