diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index fe09e6557..d04c08272 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -177,6 +177,8 @@ add_library(libslic3r STATIC PrintRegion.cpp PNGReadWrite.hpp PNGReadWrite.cpp + QuadricEdgeCollapse.cpp + QuadricEdgeCollapse.hpp Semver.cpp ShortestPath.cpp ShortestPath.hpp diff --git a/src/libslic3r/QuadricEdgeCollapse.cpp b/src/libslic3r/QuadricEdgeCollapse.cpp new file mode 100644 index 000000000..1974f5101 --- /dev/null +++ b/src/libslic3r/QuadricEdgeCollapse.cpp @@ -0,0 +1,620 @@ +#include "QuadricEdgeCollapse.hpp" +#include +#include "MutablePriorityQueue.hpp" +#include "SimplifyMeshImpl.hpp" + +using namespace Slic3r; + +// only private namespace not neccessary be in hpp +namespace QuadricEdgeCollapse { + using Vertices = std::vector; + using Triangle = stl_triangle_vertex_indices; + using Indices = std::vector; + using SymMat = SimplifyMesh::implementation::SymetricMatrix; + + struct Error + { + float value; + // range(0 .. 2), + unsigned char min_index; + Error(float value, unsigned char min_index): value(value), min_index(min_index) { + assert(min_index < 3); + } + Error() = default; + }; + using Errors = std::vector; + // merge information together - faster access during processing + struct TriangleInfo + { + Vec3f n; // normalized normal - speed up calcualtion of q and check flip + Error e; // smallest error caused by edges, identify smallest edge in triangle + TriangleInfo() = default; + bool is_deleted() const { return e.min_index > 2; } + void set_deleted() { e.min_index = 3; } + }; + using TriangleInfos = std::vector; + struct VertexInfo + { + SymMat q; // sum quadric of surround triangles + size_t start = 0, count = 0; // vertex neighbor triangles + VertexInfo() = default; + bool is_deleted() const { return count == 0; } + }; + using VertexInfos = std::vector; + struct EdgeInfo + { + size_t t_index=0; // triangle index + unsigned char edge = 0; // 0 or 1 or 2 + EdgeInfo() = default; + }; + using EdgeInfos = std::vector; + + Vec3f create_normal(const Triangle &triangle, const Vertices &vertices); + double calculate_error(size_t id_v1, size_t id_v2, SymMat & q, const Vertices &vertices); + Vec3f calculate_vertex(size_t id_v1, size_t id_v2, SymMat & q, const Vertices &vertices); + // calculate error for vertex and quadrics, triangle quadrics and triangle vertex give zero, only pozitive number + double vertex_error(const SymMat &q, const Vec3d &vertex); + SymMat create_quadric(const Triangle &t, const TriangleInfo &t_info, const Vertices &vertices); + std::tuple init(const indexed_triangle_set &its); + size_t find_triangle_index1(size_t vi, const VertexInfo& v_info, size_t ti, const EdgeInfos& e_infos, const Indices& indices); + bool is_flipped(const Vec3f &vn, const Vec3f &v1, const Vec3f &v2, const Vec3f &normal); + bool is_flipped(Vec3f &new_vertex, size_t ti0, size_t ti1, const VertexInfo& v_info, + const TriangleInfos &t_infos, const EdgeInfos &e_infos, const indexed_triangle_set &its); + + // find edge with smallest error in triangle + Error calculate_error(const Triangle& t,const Vertices &vertices, const VertexInfos& v_infos); + // subtract quadric of one triangle from triangle vertex + void sub_quadric(const Triangle &t, const TriangleInfo &t_info, VertexInfos &v_infos, const Vertices &vertices); + + void remove_triangle(EdgeInfos &e_infos, VertexInfo &v_info, size_t ti); + void change_neighbors(EdgeInfos &e_infos, VertexInfos &v_infos, size_t ti0, size_t ti1, + size_t vi0, size_t vi1, size_t vi_top0, const Triangle &t1); + + void compact(const VertexInfos &v_infos, const TriangleInfos &t_infos, const EdgeInfos &e_infos, indexed_triangle_set &its); + } + +using namespace QuadricEdgeCollapse; + +bool check_neighbors(TriangleInfos &t_infos, + Indices& indices, + VertexInfos & v_infos) +{ + std::vector t_counts(v_infos.size(), 0); + for (size_t i = 0; i < indices.size(); i++) { + TriangleInfo &t_info = t_infos[i]; + if (t_info.is_deleted()) continue; + Triangle &t = indices[i]; + for (size_t vidx : t) ++t_counts[vidx]; + } + + size_t prev_end = 0; + for (size_t i = 0; i < v_infos.size(); i++) { + VertexInfo &v_info = v_infos[i]; + if (v_info.is_deleted()) continue; + if (v_info.count != t_counts[i]) { + // no correct count + return false; + } + if (prev_end > v_info.start) { + // overlap of start + return false; + } + prev_end = v_info.start + v_info.count; + } + return true; +} + +bool check_new_vertex(const Vec3f& nv, const Vec3f& v0, const Vec3f& v1) { + float epsilon = 1.f; + for (size_t i = 0; i < 3; i++) { + if (nv[i] > (v0[i] + epsilon) && nv[i] > (v1[i] + epsilon) || + nv[i] < (v0[i] - epsilon) && nv[i] < (v1[i] - epsilon)) { + return false; + } + } + return true; +} + +bool Slic3r::its_quadric_edge_collapse(indexed_triangle_set &its, + size_t triangle_count) +{ + TriangleInfos t_infos; + VertexInfos v_infos; + EdgeInfos e_infos; + std::tie(t_infos, v_infos, e_infos) = init(its); + + static constexpr double max_error = std::numeric_limits::max(); + + auto cmp = [&t_infos](size_t vi0, size_t vi1) -> bool { + const Error &e0 = t_infos[vi0].e; + const Error &e1 = t_infos[vi1].e; + return e0.value < e1.value; + }; + // convert triangle index to priority queue index + std::vector i_convert(its.indices.size(), {0}); + auto setter = [&i_convert](size_t it, size_t index) { i_convert[it] = index; }; + MutablePriorityQueue mpq(std::move(setter), std::move(cmp)); + mpq.reserve(its.indices.size()); + for (size_t i = 0; i < its.indices.size(); i++) mpq.push(i); + + size_t actual_triangle_count = its.indices.size(); + while (actual_triangle_count > triangle_count && !mpq.empty()) { + // triangle index 0 + size_t ti0 = mpq.top(); + mpq.pop(); + TriangleInfo &t_info0 = t_infos[ti0]; + if (t_info0.is_deleted()) continue; + + Error &e = t_info0.e; + + const Triangle &t0 = its.indices[ti0]; + size_t vi0 = t0[e.min_index]; + size_t vi1 = t0[(e.min_index+1) %3]; + // Need by move of neighbor edge infos in function: change_neighbors + if (vi0 > vi1) std::swap(vi0, vi1); + VertexInfo &v_info0 = v_infos[vi0]; + VertexInfo &v_info1 = v_infos[vi1]; + assert(!v_info0.is_deleted() && !v_info1.is_deleted()); + + // new vertex position + SymMat q(v_info0.q); + q += v_info1.q; + Vec3f new_vertex0 = calculate_vertex(vi0, vi1, q, its.vertices); + assert(check_new_vertex(new_vertex0, its.vertices[vi0], its.vertices[vi1])); + // set of triangle indices that change quadric + size_t ti1 = (v_info0.count < v_info1.count)? + find_triangle_index1(vi1, v_info0, ti0, e_infos, its.indices) : + find_triangle_index1(vi0, v_info1, ti0, e_infos, its.indices) ; + + if (is_flipped(new_vertex0, ti0, ti1, v_info0, t_infos, e_infos, its) || + is_flipped(new_vertex0, ti0, ti1, v_info1, t_infos, e_infos, its)) { + // IMPROVE1: what about other edges in triangle? + // IMPROVE2: check mpq top if it is ti1 with same edge + e.value = max_error; + mpq.push(ti0); + continue; + } + + std::vector changed_triangle_indices; + changed_triangle_indices.reserve(v_info0.count + v_info1.count - 4); + + sub_quadric(t0, t_info0, v_infos, its.vertices); + TriangleInfo &t_info1 = t_infos[ti1]; + const Triangle &t1 = its.indices[ti1]; + sub_quadric(t1, t_info1, v_infos, its.vertices); + // for each vertex0 triangles + size_t v_info0_end = v_info0.start + v_info0.count; + for (size_t di = v_info0.start; di < v_info0_end; ++di) { + assert(di < e_infos.size()); + EdgeInfo &e_info = e_infos[di]; + size_t ti = e_info.t_index; + if (ti == ti0) continue; // ti0 will be deleted + if (ti == ti1) continue; // ti1 will be deleted + sub_quadric(its.indices[ti], t_infos[ti], v_infos, its.vertices); + changed_triangle_indices.emplace_back(ti); + } + + // for each vertex1 triangles + size_t v_info1_end = v_info1.start + v_info1.count; + for (size_t di = v_info1.start; di < v_info1_end; ++di) { + assert(di < e_infos.size()); + EdgeInfo &e_info = e_infos[di]; + size_t ti = e_info.t_index; + if (ti == ti0) continue; // ti0 will be deleted + if (ti == ti1) continue; // ti1 will be deleted + Triangle &t = its.indices[ti]; + sub_quadric(t, t_infos[ti], v_infos, its.vertices); + t[e_info.edge] = vi0; // change index + changed_triangle_indices.emplace_back(ti); + } + + // fix neighbors + + // vertex index of triangle 0 which is not vi0 nor vi1 + size_t vi_top0 = t0[(e.min_index + 2) % 3]; + change_neighbors(e_infos, v_infos, ti0, ti1, vi0, vi1, vi_top0, t1); + + // Change vertex + // Has to be set after subtract quadric + its.vertices[vi0] = new_vertex0; + + // add new quadrics + v_info0.q = SymMat(); // zero value + for (size_t ti : changed_triangle_indices) { + const Triangle& t = its.indices[ti]; + TriangleInfo &t_info = t_infos[ti]; + t_info.n = create_normal(t, its.vertices); // new normal + SymMat q = create_quadric(t, t_info, its.vertices); + for (const size_t vi: t) v_infos[vi].q += q; + } + + // fix errors - must be after calculate all quadric + t_info1.e.value = max_error; // not neccessary when check deleted triangles at begining + //mpq.remove(i_convert[ti1]); + for (size_t ti : changed_triangle_indices) { + const Triangle &t = its.indices[ti]; + t_infos[ti].e = calculate_error(t, its.vertices, v_infos); + mpq.update(i_convert[ti]); + } + + // set triangle(0 + 1) indices as deleted + t_info0.set_deleted(); + t_info1.set_deleted(); + // triangle counter decrementation + actual_triangle_count-=2; + + //assert(check_neighbors(t_infos, its.indices, v_infos)); + } + + // compact triangle + compact(v_infos, t_infos, e_infos, its); + return true; +} + +Vec3f QuadricEdgeCollapse::create_normal(const Triangle & triangle, + const Vertices &vertices) +{ + const Vec3f &v0 = vertices[triangle[0]]; + const Vec3f &v1 = vertices[triangle[1]]; + const Vec3f &v2 = vertices[triangle[2]]; + // n = triangle normal + Vec3f n = (v1 - v0).cross(v2 - v0); + n.normalize(); + return n; +} + +double QuadricEdgeCollapse::calculate_error(size_t id_v1, + size_t id_v2, + SymMat & q, + const Vertices &vertices) +{ + double det = q.det(0, 1, 2, 1, 4, 5, 2, 5, 7); + if (abs(det) < std::numeric_limits::epsilon()) { + // can't divide by zero + const Vec3f &v0 = vertices[id_v1]; + const Vec3f &v1 = vertices[id_v2]; + Vec3d verts[3] = {v0.cast(), v1.cast(), Vec3d()}; + verts[2] = (verts[0] + verts[1]) / 2; + double errors[] = {vertex_error(q, verts[0]), + vertex_error(q, verts[1]), + vertex_error(q, verts[2])}; + return *std::min_element(std::begin(errors), std::end(errors)); + } + + double det_1 = -1 / det; + double det_x = q.det(1, 2, 3, 4, 5, 6, 5, 7, 8); // vx = A41/det(q_delta) + double det_y = q.det(0, 2, 3, 1, 5, 6, 2, 7, 8); // vy = A42/det(q_delta) + double det_z = q.det(0, 1, 3, 1, 4, 6, 2, 5, 8); // vz = A43/det(q_delta) + Vec3d vertex(det_1 * det_x, -det_1 * det_y, det_1 * det_z); + return vertex_error(q, vertex); +} + +// similar as calculate error but focus on new vertex without calculation of error +Vec3f QuadricEdgeCollapse::calculate_vertex(size_t id_v1, + size_t id_v2, + SymMat & q, + const Vertices &vertices) +{ + double det = q.det(0, 1, 2, 1, 4, 5, 2, 5, 7); + if (abs(det) < std::numeric_limits::epsilon()) { + // can't divide by zero + const Vec3f &v0 = vertices[id_v1]; + const Vec3f &v1 = vertices[id_v2]; + Vec3d verts[3] = {v0.cast(), v1.cast(), Vec3d()}; + verts[2] = (verts[0] + verts[1]) / 2; + double errors[] = {vertex_error(q, verts[0]), + vertex_error(q, verts[1]), + vertex_error(q, verts[2])}; + auto mit = std::min_element(std::begin(errors), std::end(errors)); + return verts[mit - std::begin(errors)].cast(); + } + + double det_1 = -1 / det; + double det_x = q.det(1, 2, 3, 4, 5, 6, 5, 7, 8); // vx = A41/det(q_delta) + double det_y = q.det(0, 2, 3, 1, 5, 6, 2, 7, 8); // vy = A42/det(q_delta) + double det_z = q.det(0, 1, 3, 1, 4, 6, 2, 5, 8); // vz = A43/det(q_delta) + return Vec3f(det_1 * det_x, -det_1 * det_y, det_1 * det_z); +} + +double QuadricEdgeCollapse::vertex_error(const SymMat &q, const Vec3d &vertex) +{ + const double &x = vertex.x(), &y = vertex.y(), &z = vertex.z(); + return q[0] * x * x + 2 * q[1] * x * y + 2 * q[2] * x * z + + 2 * q[3] * x + q[4] * y * y + 2 * q[5] * y * z + + 2 * q[6] * y + q[7] * z * z + 2 * q[8] * z + q[9]; +} + +SymMat QuadricEdgeCollapse::create_quadric(const Triangle & t, + const TriangleInfo &t_info, + const Vertices & vertices) +{ + Vec3d n = t_info.n.cast(); + Vec3d v0 = vertices[t[0]].cast(); + return SymMat(n.x(), n.y(), n.z(), -n.dot(v0)); +} + +std::tuple QuadricEdgeCollapse::init( + const indexed_triangle_set &its) +{ + TriangleInfos t_infos(its.indices.size()); + VertexInfos v_infos(its.vertices.size()); + EdgeInfos e_infos(its.indices.size() * 3); + + for (size_t i = 0; i < its.indices.size(); i++) { + const Triangle &t = its.indices[i]; + TriangleInfo &t_info = t_infos[i]; + t_info.n = create_normal(t, its.vertices); + SymMat q = create_quadric(t, t_info, its.vertices); + for (size_t e = 0; e < 3; e++) { + VertexInfo &v_info = v_infos[t[e]]; + v_info.q += q; + ++v_info.count; // triangle count + } + } + + // set offseted starts + size_t triangle_start = 0; + for (VertexInfo &v_info : v_infos) { + v_info.start = triangle_start; + triangle_start += v_info.count; + // set filled vertex to zero + v_info.count = 0; + } + assert(its.indices.size() * 3 == triangle_start); + + // calc error + create reference + for (size_t i = 0; i < its.indices.size(); i++) { + const Triangle &t = its.indices[i]; + TriangleInfo &t_info = t_infos[i]; + t_info.e = calculate_error(t, its.vertices, v_infos); + for (size_t j = 0; j < 3; ++j) { + VertexInfo &v_info = v_infos[t[j]]; + size_t ei = v_info.start + v_info.count; + assert(ei < e_infos.size()); + EdgeInfo &e_info = e_infos[ei]; + e_info.t_index = i; + e_info.edge = j; + ++v_info.count; + } + } + return {t_infos, v_infos, e_infos}; +} + +size_t QuadricEdgeCollapse::find_triangle_index1(size_t vi, + const VertexInfo &v_info, + size_t ti0, + const EdgeInfos & e_infos, + const Indices &indices) +{ + size_t end = v_info.start + v_info.count; + for (size_t ei = v_info.start; ei < end; ++ei) { + const EdgeInfo &e_info = e_infos[ei]; + if (e_info.t_index == ti0) continue; + const Triangle& t = indices[e_info.t_index]; + if (t[(e_info.edge + 1) % 3] == vi || + t[(e_info.edge + 2) % 3] == vi) + return e_info.t_index; + } + // triangle0 is on border and do NOT have twin edge + assert(false); + return -1; +} + + +bool QuadricEdgeCollapse::is_flipped(const Vec3f &vn, + const Vec3f &v1, + const Vec3f &v2, + const Vec3f &normal) +{ + static const float thr_pos = 1.0f - std::numeric_limits::epsilon(); + static const float thr_neg = -thr_pos; + static const float dot_thr = 0.2f; // Value from simplify mesh + + Vec3f d1 = v1 - vn; + d1.normalize(); + Vec3f d2 = v2 - vn; + d2.normalize(); + + float dot = d1.dot(d2); + if (dot > thr_pos || dot < thr_neg) return true; + + // IMPROVE: propagate new normal + Vec3f n = d1.cross(d2); + n.normalize(); + return n.dot(normal) < dot_thr; +} + + +bool QuadricEdgeCollapse::is_flipped(Vec3f & new_vertex, + size_t ti0, + size_t ti1, + const VertexInfo & v_info, + const TriangleInfos &t_infos, + const EdgeInfos & e_infos, + const indexed_triangle_set &its) +{ + // for each vertex triangles + size_t v_info_end = v_info.start + v_info.count; + for (size_t ei = v_info.start; ei < v_info_end; ++ei) { + assert(ei < e_infos.size()); + const EdgeInfo &e_info = e_infos[ei]; + if (e_info.t_index == ti0) continue; // ti0 will be deleted + if (e_info.t_index == ti1) continue; // ti1 will be deleted + const Triangle &t = its.indices[e_info.t_index]; + const Vec3f &normal = t_infos[e_info.t_index].n; + const Vec3f &vf = its.vertices[t[(e_info.edge + 1) % 3]]; + const Vec3f &vs = its.vertices[t[(e_info.edge + 2) % 3]]; + if (is_flipped(new_vertex, vf, vs, normal)) + return true; + } + return false; +} + +Error QuadricEdgeCollapse::calculate_error(const Triangle & t, + const Vertices & vertices, + const VertexInfos &v_infos) +{ + Vec3d error; + for (size_t j = 0; j < 3; ++j) { + size_t j2 = (j == 2) ? 0 : (j + 1); + size_t vi0 = t[j]; + size_t vi1 = t[j2]; + SymMat q(v_infos[vi0].q); // copy + q += v_infos[vi1].q; + error[j] = calculate_error(vi0, vi1, q, vertices); + } + unsigned char min_index = (error[0] < error[1]) ? + ((error[0] < error[2]) ? 0 : 2) : + ((error[1] < error[2]) ? 1 : 2); + return Error(static_cast(error[min_index]), min_index); +} + +void QuadricEdgeCollapse::sub_quadric(const Triangle &t, + const TriangleInfo & t_info, + VertexInfos &v_infos, + const Vertices &vertices) +{ + SymMat quadric = create_quadric(t, t_info, vertices); + for (size_t vi: t) v_infos[vi].q -= quadric; +} + +void QuadricEdgeCollapse::remove_triangle(EdgeInfos & e_infos, + VertexInfo &v_info, + size_t ti) +{ + auto e_info = e_infos.begin() + v_info.start; + auto e_info_end = e_info + v_info.count - 1; + for (; e_info != e_info_end; ++e_info) { + if (e_info->t_index == ti) { + *e_info = *e_info_end; + --v_info.count; + return; + } + } + assert(e_info_end->t_index == ti); + // last triangle is ti + --v_info.count; +} + +void QuadricEdgeCollapse::change_neighbors(EdgeInfos & e_infos, + VertexInfos & v_infos, + size_t ti0, + size_t ti1, + size_t vi0, + size_t vi1, + size_t vi_top0, + const Triangle &t1) +{ + // have to copy Edge info from higher vertex index into smaller + assert(vi0 < vi1); + + // vertex index of triangle 1 which is not vi0 nor vi1 + size_t vi_top1 = t1[0]; + if (vi_top1 == vi0 || vi_top1 == vi1) { + vi_top1 = (t1[1] == vi0 || t1[1] == vi1)? t1[2] : t1[1]; + } + + remove_triangle(e_infos, v_infos[vi_top0], ti0); + remove_triangle(e_infos, v_infos[vi_top1], ti1); + + VertexInfo &v_info0 = v_infos[vi0]; + VertexInfo &v_info1 = v_infos[vi1]; + + size_t new_triangle_count = v_info0.count + v_info1.count - 4; + remove_triangle(e_infos, v_info0, ti0); + remove_triangle(e_infos, v_info0, ti1); + + // copy second's edge infos out of e_infos, to free size + EdgeInfos e_infos1; + e_infos1.reserve(v_info1.count - 2); + size_t v_info_s_end = v_info1.start + v_info1.count; + for (size_t ei = v_info1.start; ei < v_info_s_end; ++ei) { + const EdgeInfo &e_info = e_infos[ei]; + if (e_info.t_index == ti0) continue; + if (e_info.t_index == ti1) continue; + e_infos1.emplace_back(e_info); + } + v_info1.count = 0; + + size_t need = (new_triangle_count < v_info0.count)? 0: + (new_triangle_count - v_info0.count); + + size_t act_vi = vi0 + 1; + VertexInfo *act_v_info = &v_infos[act_vi]; + size_t act_start = act_v_info->start; + size_t last_end = v_info0.start + v_info0.count; + + struct CopyEdgeInfo + { + size_t start; + size_t count; + unsigned char move; + CopyEdgeInfo(size_t start, size_t count, unsigned char move) + : start(start), count(count), move(move) + {} + }; + std::vector c_infos; + c_infos.reserve(need); + while (true) { + size_t save = act_start - last_end; + if (save > 0) { + if (save >= need) break; + need -= save; + c_infos.emplace_back(act_v_info->start, act_v_info->count, need); + } else { + c_infos.back().count += act_v_info->count; + } + last_end = act_v_info->start + act_v_info->count; + act_v_info->start += need; + ++act_vi; + if (act_vi < v_infos.size()) { + act_v_info = &v_infos[act_vi]; + act_start = act_v_info->start; + } else + act_start = e_infos.size(); // fix for edge between last two triangles + } + + // copy by c_infos + for (size_t i = c_infos.size(); i > 0; --i) { + const CopyEdgeInfo &c_info = c_infos[i - 1]; + for (size_t ei = c_info.start + c_info.count - 1; ei >= c_info.start; --ei) + e_infos[ei + c_info.move] = e_infos[ei]; // copy + } + + // copy triangle from first info into second + for (size_t ei_s = 0; ei_s < e_infos1.size(); ++ei_s) { + size_t ei_f = v_info0.start + v_info0.count; + e_infos[ei_f] = e_infos1[ei_s]; // copy + ++v_info0.count; + } +} + +void QuadricEdgeCollapse::compact(const VertexInfos & v_infos, + const TriangleInfos & t_infos, + const EdgeInfos & e_infos, + indexed_triangle_set &its) +{ + size_t vi_new = 0; + for (size_t vi = 0; vi < v_infos.size(); vi++) { + const VertexInfo &v_info = v_infos[vi]; + if (v_info.is_deleted()) continue; // deleted + size_t e_info_end = v_info.start + v_info.count; + for (size_t ei = v_info.start; ei < e_info_end; ei++) { + const EdgeInfo &e_info = e_infos[ei]; + // change vertex index + its.indices[e_info.t_index][e_info.edge] = vi_new; + } + // compact vertices + its.vertices[vi_new++] = its.vertices[vi]; + } + // remove vertices tail + its.vertices.erase(its.vertices.begin() + vi_new, its.vertices.end()); + + size_t ti_new = 0; + for (size_t ti = 0; ti < t_infos.size(); ti++) { + const TriangleInfo &t_info = t_infos[ti]; + if (t_info.is_deleted()) continue; + its.indices[ti_new++] = its.indices[ti]; + } + its.indices.erase(its.indices.begin() + ti_new, its.indices.end()); +} \ No newline at end of file diff --git a/src/libslic3r/QuadricEdgeCollapse.hpp b/src/libslic3r/QuadricEdgeCollapse.hpp new file mode 100644 index 000000000..d435656b2 --- /dev/null +++ b/src/libslic3r/QuadricEdgeCollapse.hpp @@ -0,0 +1,17 @@ +// paper: https://people.eecs.berkeley.edu/~jrs/meshpapers/GarlandHeckbert2.pdf +// sum up: https://users.csc.calpoly.edu/~zwood/teaching/csc570/final06/jseeba/ +// inspiration: https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification + +#include "TriangleMesh.hpp" + +namespace Slic3r { + +/// +/// Simplify mesh by Quadric metric +/// +/// IN/OUT triangle mesh to be simplified. +/// wanted triangle count. +/// TRUE on success otherwise FALSE +bool its_quadric_edge_collapse(indexed_triangle_set &its, size_t triangle_count); + +} // namespace Slic3r diff --git a/src/libslic3r/SimplifyMeshImpl.hpp b/src/libslic3r/SimplifyMeshImpl.hpp index 4b6b0f5cb..0b3b9a2f6 100644 --- a/src/libslic3r/SimplifyMeshImpl.hpp +++ b/src/libslic3r/SimplifyMeshImpl.hpp @@ -107,7 +107,7 @@ public: // Determinant T det(int a11, int a12, int a13, int a21, int a22, int a23, - int a31, int a32, int a33) + int a31, int a32, int a33) const { T det = m[a11] * m[a22] * m[a33] + m[a13] * m[a21] * m[a32] + m[a12] * m[a23] * m[a31] - m[a13] * m[a22] * m[a31] - @@ -121,7 +121,13 @@ public: for (size_t i = 0; i < N; ++i) m[i] += n[i]; return *this; } - + + const SymetricMatrix &operator-=(const SymetricMatrix &n) + { + for (size_t i = 0; i < N; ++i) m[i] -= n[i]; + return *this; + } + SymetricMatrix operator+(const SymetricMatrix& n) { SymetricMatrix self = *this; diff --git a/tests/libslic3r/test_indexed_triangle_set.cpp b/tests/libslic3r/test_indexed_triangle_set.cpp index 4c1128345..0ec7e3bd0 100644 --- a/tests/libslic3r/test_indexed_triangle_set.cpp +++ b/tests/libslic3r/test_indexed_triangle_set.cpp @@ -4,6 +4,8 @@ #include "libslic3r/TriangleMesh.hpp" +using namespace Slic3r; + TEST_CASE("Split empty mesh", "[its_split][its]") { using namespace Slic3r; @@ -100,3 +102,41 @@ TEST_CASE("Split two watertight meshes", "[its_split][its]") { debug_write_obj(res, "parts_watertight"); } +#include +TEST_CASE("Reduce one edge by Quadric Edge Collapse", "[its]") +{ + indexed_triangle_set its; + its.vertices = {Vec3f(-1.f, 0.f, 0.f), Vec3f(0.f, 1.f, 0.f), + Vec3f(1.f, 0.f, 0.f), Vec3f(0.f, 0.f, 1.f), + // vertex to be removed + Vec3f(0.9f, .1f, -.1f)}; + its.indices = {Vec3i(1, 0, 3), Vec3i(2, 1, 3), Vec3i(0, 2, 3), + Vec3i(0, 1, 4), Vec3i(1, 2, 4), Vec3i(2, 0, 4)}; + // edge to remove is between vertices 2 and 4 on trinagles 4 and 5 + + indexed_triangle_set its_ = its; // copy + // its_write_obj(its, "tetrhedron_in.obj"); + size_t wanted_count = its.indices.size() - 1; + CHECK(its_quadric_edge_collapse(its, wanted_count)); + // its_write_obj(its, "tetrhedron_out.obj"); + CHECK(its.indices.size() == 4); + CHECK(its.vertices.size() == 4); + + for (size_t i = 0; i < 3; i++) { + CHECK(its.indices[i] == its_.indices[i]); + } + + for (size_t i = 0; i < 4; i++) { + if (i == 2) continue; + CHECK(its.vertices[i] == its_.vertices[i]); + } + + const Vec3f &v = its.vertices[2]; // new vertex + const Vec3f &v2 = its_.vertices[2]; // moved vertex + const Vec3f &v4 = its_.vertices[4]; // removed vertex + for (size_t i = 0; i < 3; i++) { + bool is_between = (v[i] < v4[i] && v[i] > v2[i]) || + (v[i] > v4[i] && v[i] < v2[i]); + CHECK(is_between); + } +} \ No newline at end of file