From acc3d872b3164614fa7cb98b3dd7c5523796fc4d Mon Sep 17 00:00:00 2001 From: pavel Date: Thu, 25 Jun 2020 01:55:04 +0200 Subject: [PATCH] added quadric and memoryless simplification --- CMakeLists.txt | 2 + Graph.hpp | 32 +-- Math.hpp | 7 +- Matrix.hpp | 175 +++++++++++----- MemorylessDecimator.hpp | 430 ++++++++++++++++++++++++++++++++++++++++ PlyReader.hpp | 11 +- PlyWriter.hpp | 4 +- QuadricDecimator.hpp | 79 ++++++++ Simplification.hpp | 137 +++++++------ TriangleMesh.hpp | 22 ++ 10 files changed, 770 insertions(+), 129 deletions(-) create mode 100644 MemorylessDecimator.hpp create mode 100644 QuadricDecimator.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f95739d..afd4110 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,6 +32,8 @@ add_executable(pvl TriangleMesh.hpp Refinement.hpp Simplification.hpp + QuadricDecimator.hpp + MemorylessDecimator.hpp Range.hpp Graph.hpp Graph.inl.hpp UniformGrid.hpp diff --git a/Graph.hpp b/Graph.hpp index 904896a..af5f7cd 100644 --- a/Graph.hpp +++ b/Graph.hpp @@ -655,18 +655,22 @@ class Graph { } struct CollapseContext { - HalfEdgeHandle edge; // halfedge from "remaining" to "removed" + HalfEdgeHandle edge; VertexHandle removed; VertexHandle remaining; FaceHandle left; // left of the halfedge FaceHandle right; // right of the halfedge - CollapseContext(const Graph& graph, const HalfEdgeHandle heh) - : edge(heh) { + CollapseContext(const Graph& graph, const EdgeHandle eh) + : edge(graph.halfEdge(eh)) { removed = graph.to(edge); remaining = graph.from(edge); left = graph.left(edge); - right = graph.right(edge); + if (!graph.boundary(edge)) { + right = graph.right(edge); + } else { + right = FaceHandle(-1); + } } }; @@ -675,7 +679,7 @@ class Graph { void collapse(HalfEdgeHandle heh) { PVL_ASSERT(collapseAllowed(heh)); // check(); - CollapseContext context(*this, heh); + CollapseContext context(*this, edge(heh)); HalfEdgeHandle oheh = opposite(heh); /*if (vertices_[context.remaining].edge == heh) { @@ -762,12 +766,14 @@ class Graph { std::vector is; std::set_intersection( ring1.begin(), ring1.end(), ring2.begin(), ring2.end(), std::back_inserter(is)); - if (is.size() != 2) { - return false; - } - VertexHandle vl = to(next(eh)); - VertexHandle vr = to(next(opposite(eh))); - return ((is[0] == vl && is[1] == vr) || (is[1] == vl && is[0] == vr)); + PVL_ASSERT(is.size() >= 2); + return is.size() == 2; + // if (is.size() != 2) { + // return false; + //} + // VertexHandle vl = to(next(eh)); + // VertexHandle vr = to(next(opposite(eh))); + // return ((is[0] == vl && is[1] == vr) || (is[1] == vl && is[0] == vr)); } void remove(FaceHandle fh) { @@ -797,8 +803,8 @@ class Graph { return halfEdges_[heh].left == -1; } - bool removed(EdgeHandle heh) { - return halfEdges_[heh].left == -1; + bool removed(EdgeHandle eh) { + return removed(halfEdge(eh)); } void collectGarbage() {} diff --git a/Math.hpp b/Math.hpp index 382a813..bba2f6f 100644 --- a/Math.hpp +++ b/Math.hpp @@ -31,13 +31,18 @@ struct Pow { return value * value; } }; - template struct Pow { T operator()(const T& value) { return value * value * value; } }; +template +struct Pow { + T operator()(const T& value) { + return sqr(value * value); + } +}; template T pow(const T& value) { diff --git a/Matrix.hpp b/Matrix.hpp index 90d04f1..3b1e73d 100644 --- a/Matrix.hpp +++ b/Matrix.hpp @@ -32,12 +32,16 @@ class Matrix { return rows_[r][c]; } - const Vector& row(const int idx) const { PVL_ASSERT(idx < Rows); return rows_[idx]; } + Vector& row(const int idx) { + PVL_ASSERT(idx < Rows); + return rows_[idx]; + } + Vector column(const int idx) const { PVL_ASSERT(idx < Cols); Vector col; @@ -89,6 +93,27 @@ class Matrix { return m; } + Matrix operator*(const T f) const { + Matrix m; + for (int j = 0; j < Rows; ++j) { + for (int i = 0; i < Cols; ++i) { + m(i, j) = (*this)(i, j) * f; + } + } + return m; + } + + Matrix operator/(const T f) const { + PVL_ASSERT(f != 0); + Matrix m; + for (int j = 0; j < Rows; ++j) { + for (int i = 0; i < Cols; ++i) { + m(i, j) = (*this)(i, j) / f; + } + } + return m; + } + Matrix& operator+=(const Matrix& other) { *this = *this + other; return *this; @@ -114,16 +139,26 @@ class Matrix { using Mat33f = Matrix; using Mat44f = Matrix; -inline Mat33f outerProd(const Vec3f& v1, const Vec3f& v2) { - Mat33f res; - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { +template +Matrix outerProd(const Vector& v1, const Vector& v2) { + Matrix res; + for (int i = 0; i < N; ++i) { + for (int j = 0; j < M; ++j) { res(i, j) = v1[i] * v2[j]; } } return res; } +template +inline T trace(const Matrix& m) { + T tr = 0.; + for (int i = 0; i < N; ++i) { + tr += m(i, i); + } + return tr; +} + inline Mat33f invert(const Mat33f& m) { double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) + @@ -144,57 +179,97 @@ inline Mat33f invert(const Mat33f& m) { return minv; } -inline Mat44f invert(const Mat44f& m) { - Mat44f minv; - minv(0, 0) = m(1, 1) * m(2, 2) * m(3, 3) - m(1, 1) * m(3, 2) * m(2, 3) - m(1, 2) * m(2, 1) * m(3, 3) + - m(1, 2) * m(3, 1) * m(2, 3) + m(1, 3) * m(2, 1) * m(3, 2) - m(1, 3) * m(3, 1) * m(2, 2); - - minv(0, 1) = -m(0, 1) * m(2, 2) * m(3, 3) + m(0, 1) * m(3, 2) * m(2, 3) + m(0, 2) * m(2, 1) * m(3, 3) - - m(0, 2) * m(3, 1) * m(2, 3) - m(0, 3) * m(2, 1) * m(3, 2) + m(0, 3) * m(3, 1) * m(2, 2); - - minv(0, 2) = m(0, 1) * m(1, 2) * m(3, 3) - m(0, 1) * m(3, 2) * m(1, 3) - m(0, 2) * m(1, 1) * m(3, 3) + - m(0, 2) * m(3, 1) * m(1, 3) + m(0, 3) * m(1, 1) * m(3, 2) - m(0, 3) * m(3, 1) * m(1, 2); +inline float determinant(const Mat44f& m) { + Mat44f minv; /// \todo + minv(0, 0) = m(1, 1) * m(2, 2) * m(3, 3) - m(1, 1) * m(3, 2) * m(2, 3) - + m(1, 2) * m(2, 1) * m(3, 3) + m(1, 2) * m(3, 1) * m(2, 3) + + m(1, 3) * m(2, 1) * m(3, 2) - m(1, 3) * m(3, 1) * m(2, 2); - minv(0, 3) = -m(0, 1) * m(1, 2) * m(2, 3) + m(0, 1) * m(2, 2) * m(1, 3) + m(0, 2) * m(1, 1) * m(2, 3) - - m(0, 2) * m(2, 1) * m(1, 3) - m(0, 3) * m(1, 1) * m(2, 2) + m(0, 3) * m(2, 1) * m(1, 2); + minv(0, 1) = -m(0, 1) * m(2, 2) * m(3, 3) + m(0, 1) * m(3, 2) * m(2, 3) + + m(0, 2) * m(2, 1) * m(3, 3) - m(0, 2) * m(3, 1) * m(2, 3) - + m(0, 3) * m(2, 1) * m(3, 2) + m(0, 3) * m(3, 1) * m(2, 2); - minv(1, 0) = -m(1, 0) * m(2, 2) * m(3, 3) + m(1, 0) * m(3, 2) * m(2, 3) + m(1, 2) * m(2, 0) * m(3, 3) - - m(1, 2) * m(3, 0) * m(2, 3) - m(1, 3) * m(2, 0) * m(3, 2) + m(1, 3) * m(3, 0) * m(2, 2); + minv(0, 2) = m(0, 1) * m(1, 2) * m(3, 3) - m(0, 1) * m(3, 2) * m(1, 3) - + m(0, 2) * m(1, 1) * m(3, 3) + m(0, 2) * m(3, 1) * m(1, 3) + + m(0, 3) * m(1, 1) * m(3, 2) - m(0, 3) * m(3, 1) * m(1, 2); - minv(1, 1) = m(0, 0) * m(2, 2) * m(3, 3) - m(0, 0) * m(3, 2) * m(2, 3) - m(0, 2) * m(2, 0) * m(3, 3) + - m(0, 2) * m(3, 0) * m(2, 3) + m(0, 3) * m(2, 0) * m(3, 2) - m(0, 3) * m(3, 0) * m(2, 2); + minv(0, 3) = -m(0, 1) * m(1, 2) * m(2, 3) + m(0, 1) * m(2, 2) * m(1, 3) + + m(0, 2) * m(1, 1) * m(2, 3) - m(0, 2) * m(2, 1) * m(1, 3) - + m(0, 3) * m(1, 1) * m(2, 2) + m(0, 3) * m(2, 1) * m(1, 2); - minv(1, 2) = -m(0, 0) * m(1, 2) * m(3, 3) + m(0, 0) * m(3, 2) * m(1, 3) + m(0, 2) * m(1, 0) * m(3, 3) - - m(0, 2) * m(3, 0) * m(1, 3) - m(0, 3) * m(1, 0) * m(3, 2) + m(0, 3) * m(3, 0) * m(1, 2); - - minv(1, 3) = m(0, 0) * m(1, 2) * m(2, 3) - m(0, 0) * m(2, 2) * m(1, 3) - m(0, 2) * m(1, 0) * m(2, 3) + - m(0, 2) * m(2, 0) * m(1, 3) + m(0, 3) * m(1, 0) * m(2, 2) - m(0, 3) * m(2, 0) * m(1, 2); - - minv(2, 0) = m(1, 0) * m(2, 1) * m(3, 3) - m(1, 0) * m(3, 1) * m(2, 3) - m(1, 1) * m(2, 0) * m(3, 3) + - m(1, 1) * m(3, 0) * m(2, 3) + m(1, 3) * m(2, 0) * m(3, 1) - m(1, 3) * m(3, 0) * m(2, 1); - - minv(2, 1) = -m(0, 0) * m(2, 1) * m(3, 3) + m(0, 0) * m(3, 1) * m(2, 3) + m(0, 1) * m(2, 0) * m(3, 3) - - m(0, 1) * m(3, 0) * m(2, 3) - m(0, 3) * m(2, 0) * m(3, 1) + m(0, 3) * m(3, 0) * m(2, 1); - - minv(2, 2) = m(0, 0) * m(1, 1) * m(3, 3) - m(0, 0) * m(3, 1) * m(1, 3) - m(0, 1) * m(1, 0) * m(3, 3) + - m(0, 1) * m(3, 0) * m(1, 3) + m(0, 3) * m(1, 0) * m(3, 1) - m(0, 3) * m(3, 0) * m(1, 1); - - minv(2, 3) = -m(0, 0) * m(1, 1) * m(2, 3) + m(0, 0) * m(2, 1) * m(1, 3) + m(0, 1) * m(1, 0) * m(2, 3) - - m(0, 1) * m(2, 0) * m(1, 3) - m(0, 3) * m(1, 0) * m(2, 1) + m(0, 3) * m(2, 0) * m(1, 1); - - minv(3, 0) = -m(1, 0) * m(2, 1) * m(3, 2) + m(1, 0) * m(3, 1) * m(2, 2) + m(1, 1) * m(2, 0) * m(3, 2) - - m(1, 1) * m(3, 0) * m(2, 2) - m(1, 2) * m(2, 0) * m(3, 1) + m(1, 2) * m(3, 0) * m(2, 1); - - minv(3, 1) = m(0, 0) * m(2, 1) * m(3, 2) - m(0, 0) * m(3, 1) * m(2, 2) - m(0, 1) * m(2, 0) * m(3, 2) + - m(0, 1) * m(3, 0) * m(2, 2) + m(0, 2) * m(2, 0) * m(3, 1) - m(0, 2) * m(3, 0) * m(2, 1); + const float det = m(0, 0) * minv(0, 0) + m(1, 0) * minv(0, 1) + m(2, 0) * minv(0, 2) + + m(3, 0) * minv(0, 3); + return det; +} - minv(3, 2) = -m(0, 0) * m(1, 1) * m(3, 2) + m(0, 0) * m(3, 1) * m(1, 2) + m(0, 1) * m(1, 0) * m(3, 2) - - m(0, 1) * m(3, 0) * m(1, 2) - m(0, 2) * m(1, 0) * m(3, 1) + m(0, 2) * m(3, 0) * m(1, 1); - minv(3, 3) = m(0, 0) * m(1, 1) * m(2, 2) - m(0, 0) * m(2, 1) * m(1, 2) - m(0, 1) * m(1, 0) * m(2, 2) + - m(0, 1) * m(2, 0) * m(1, 2) + m(0, 2) * m(1, 0) * m(2, 1) - m(0, 2) * m(2, 0) * m(1, 1); - const float det = - m(0, 0) * minv(0, 0) + m(1, 0) * minv(0, 1) + m(2, 0) * minv(0, 2) + m(3, 0) * minv(0, 3); +inline Mat44f invert(const Mat44f& m) { + Mat44f minv; + minv(0, 0) = m(1, 1) * m(2, 2) * m(3, 3) - m(1, 1) * m(3, 2) * m(2, 3) - + m(1, 2) * m(2, 1) * m(3, 3) + m(1, 2) * m(3, 1) * m(2, 3) + + m(1, 3) * m(2, 1) * m(3, 2) - m(1, 3) * m(3, 1) * m(2, 2); + + minv(0, 1) = -m(0, 1) * m(2, 2) * m(3, 3) + m(0, 1) * m(3, 2) * m(2, 3) + + m(0, 2) * m(2, 1) * m(3, 3) - m(0, 2) * m(3, 1) * m(2, 3) - + m(0, 3) * m(2, 1) * m(3, 2) + m(0, 3) * m(3, 1) * m(2, 2); + + minv(0, 2) = m(0, 1) * m(1, 2) * m(3, 3) - m(0, 1) * m(3, 2) * m(1, 3) - + m(0, 2) * m(1, 1) * m(3, 3) + m(0, 2) * m(3, 1) * m(1, 3) + + m(0, 3) * m(1, 1) * m(3, 2) - m(0, 3) * m(3, 1) * m(1, 2); + + minv(0, 3) = -m(0, 1) * m(1, 2) * m(2, 3) + m(0, 1) * m(2, 2) * m(1, 3) + + m(0, 2) * m(1, 1) * m(2, 3) - m(0, 2) * m(2, 1) * m(1, 3) - + m(0, 3) * m(1, 1) * m(2, 2) + m(0, 3) * m(2, 1) * m(1, 2); + + minv(1, 0) = -m(1, 0) * m(2, 2) * m(3, 3) + m(1, 0) * m(3, 2) * m(2, 3) + + m(1, 2) * m(2, 0) * m(3, 3) - m(1, 2) * m(3, 0) * m(2, 3) - + m(1, 3) * m(2, 0) * m(3, 2) + m(1, 3) * m(3, 0) * m(2, 2); + + minv(1, 1) = m(0, 0) * m(2, 2) * m(3, 3) - m(0, 0) * m(3, 2) * m(2, 3) - + m(0, 2) * m(2, 0) * m(3, 3) + m(0, 2) * m(3, 0) * m(2, 3) + + m(0, 3) * m(2, 0) * m(3, 2) - m(0, 3) * m(3, 0) * m(2, 2); + + minv(1, 2) = -m(0, 0) * m(1, 2) * m(3, 3) + m(0, 0) * m(3, 2) * m(1, 3) + + m(0, 2) * m(1, 0) * m(3, 3) - m(0, 2) * m(3, 0) * m(1, 3) - + m(0, 3) * m(1, 0) * m(3, 2) + m(0, 3) * m(3, 0) * m(1, 2); + + minv(1, 3) = m(0, 0) * m(1, 2) * m(2, 3) - m(0, 0) * m(2, 2) * m(1, 3) - + m(0, 2) * m(1, 0) * m(2, 3) + m(0, 2) * m(2, 0) * m(1, 3) + + m(0, 3) * m(1, 0) * m(2, 2) - m(0, 3) * m(2, 0) * m(1, 2); + + minv(2, 0) = m(1, 0) * m(2, 1) * m(3, 3) - m(1, 0) * m(3, 1) * m(2, 3) - + m(1, 1) * m(2, 0) * m(3, 3) + m(1, 1) * m(3, 0) * m(2, 3) + + m(1, 3) * m(2, 0) * m(3, 1) - m(1, 3) * m(3, 0) * m(2, 1); + + minv(2, 1) = -m(0, 0) * m(2, 1) * m(3, 3) + m(0, 0) * m(3, 1) * m(2, 3) + + m(0, 1) * m(2, 0) * m(3, 3) - m(0, 1) * m(3, 0) * m(2, 3) - + m(0, 3) * m(2, 0) * m(3, 1) + m(0, 3) * m(3, 0) * m(2, 1); + + minv(2, 2) = m(0, 0) * m(1, 1) * m(3, 3) - m(0, 0) * m(3, 1) * m(1, 3) - + m(0, 1) * m(1, 0) * m(3, 3) + m(0, 1) * m(3, 0) * m(1, 3) + + m(0, 3) * m(1, 0) * m(3, 1) - m(0, 3) * m(3, 0) * m(1, 1); + + minv(2, 3) = -m(0, 0) * m(1, 1) * m(2, 3) + m(0, 0) * m(2, 1) * m(1, 3) + + m(0, 1) * m(1, 0) * m(2, 3) - m(0, 1) * m(2, 0) * m(1, 3) - + m(0, 3) * m(1, 0) * m(2, 1) + m(0, 3) * m(2, 0) * m(1, 1); + + minv(3, 0) = -m(1, 0) * m(2, 1) * m(3, 2) + m(1, 0) * m(3, 1) * m(2, 2) + + m(1, 1) * m(2, 0) * m(3, 2) - m(1, 1) * m(3, 0) * m(2, 2) - + m(1, 2) * m(2, 0) * m(3, 1) + m(1, 2) * m(3, 0) * m(2, 1); + + minv(3, 1) = m(0, 0) * m(2, 1) * m(3, 2) - m(0, 0) * m(3, 1) * m(2, 2) - + m(0, 1) * m(2, 0) * m(3, 2) + m(0, 1) * m(3, 0) * m(2, 2) + + m(0, 2) * m(2, 0) * m(3, 1) - m(0, 2) * m(3, 0) * m(2, 1); + + minv(3, 2) = -m(0, 0) * m(1, 1) * m(3, 2) + m(0, 0) * m(3, 1) * m(1, 2) + + m(0, 1) * m(1, 0) * m(3, 2) - m(0, 1) * m(3, 0) * m(1, 2) - + m(0, 2) * m(1, 0) * m(3, 1) + m(0, 2) * m(3, 0) * m(1, 1); + + minv(3, 3) = m(0, 0) * m(1, 1) * m(2, 2) - m(0, 0) * m(2, 1) * m(1, 2) - + m(0, 1) * m(1, 0) * m(2, 2) + m(0, 1) * m(2, 0) * m(1, 2) + + m(0, 2) * m(1, 0) * m(2, 1) - m(0, 2) * m(2, 0) * m(1, 1); + const float det = m(0, 0) * minv(0, 0) + m(1, 0) * minv(0, 1) + m(2, 0) * minv(0, 2) + + m(3, 0) * minv(0, 3); const float norm = 1.f / det; for (int i = 0; i < 4; i++) { diff --git a/MemorylessDecimator.hpp b/MemorylessDecimator.hpp new file mode 100644 index 0000000..de6a07d --- /dev/null +++ b/MemorylessDecimator.hpp @@ -0,0 +1,430 @@ +#pragma once + +#include "Matrix.hpp" +#include "TriangleMesh.hpp" + +namespace Pvl { + +static constexpr double ALPHA = 1. * M_PI / 180.; // 1 degree +static double SQR_COS_ALPHA = sqr(std::cos(ALPHA)); +static double SQR_SIN_ALPHA = sqr(std::sin(ALPHA)); + +static constexpr double MAX_PHI = 90. * M_PI / 180.; +static double MIN_COS_PHI = std::cos(MAX_PHI); + +/*class OnScopeExit : boost::noncopyable { + std::function closure_; + +public: + OnScopeExit(std::function&& closure) + : closure_(std::move(closure)) {} + + ~OnScopeExit() { + closure_(); + } +};*/ + +class LindstromTurkConstraints { + Mat33f A_; + Vec3f b_; + int n_ = 0; + +public: + bool addConstraint(const Vec3f& p, const double b) { + if (full()) { + return false; + } + + if (isCompatible(p)) { + const double distSqr = dot(p, p); + const double dist = sqrt(distSqr); + const Vec3f pn = p / dist; + const double bn = b / dist; + + A_.row(n_) = pn; + b_[n_] = bn; + n_++; + return true; + } else { + std::cout << "Constraint not compatible" << std::endl; + return false; + } + } + + /// Used for problems formulated as quadratic minimization, + /// i.e. finds point where gradient Hv + c = 0. + int addQuadraticConstraint(const Mat33f& H, const Vec3f& c) { + int added = 0; + switch (n_) { + case 0: + added += int(addConstraint(H.row(0), -c[0])); + added += int(addConstraint(H.row(1), -c[1])); + added += int(addConstraint(H.row(2), -c[2])); + break; + case 1: { + const Vec3f r0 = A_.row(0); + std::array abs_r0 = { + std::abs(r0[0]), std::abs(r0[1]), std::abs(r0[2]) + }; + Vec3f q0; + const int maxIdx = + int(std::max_element(abs_r0.begin(), abs_r0.end()) - abs_r0.begin()); + switch (maxIdx) { + case 0: + q0 = { -r0[2] / r0[0], 0., 1. }; + break; + case 1: + q0 = { 0., -r0[2] / r0[1], 1. }; + break; + case 2: + q0 = { 1., 0., -r0[0] / r0[2] }; + break; + } + const Vec3f q1 = cross(r0, q0); + const Vec3f p0 = H.transform(q0); + const Vec3f p1 = H.transform(q1); + const double b0 = -dot(q0, c); + const double b1 = -dot(q1, c); + + added += int(addConstraint(p0, b0)); + added += int(addConstraint(p1, b1)); + break; + } + case 2: { + const Vec3f r0 = A_.row(0); + const Vec3f r1 = A_.row(1); + const Vec3f n = cross(r0, r1); + const Vec3f p = H.transform(n); + const double b = -dot(n, c); + added += int(addConstraint(p, b)); + break; + } + default: + std::cout << "Overconstrained linear system, skipping" << std::endl; + break; + } + return added; + } + + Vec3f getPlacement() const { + const Mat33f Ainv = invert(A_); + return Ainv.transform(b_); + } + + bool full() const { + return n_ == 3; + } + + int count() const { + return n_; + } + +private: + /// Contraint compatibility rules (see Sec. 4.2 in LT paper) + bool isCompatible(const Vec3f& p) const { + const double distSqr = normSqr(p); + if (distSqr < 1.e-20) { + return false; + } + switch (n_) { + case 0: + return true; + case 1: { + const Vec3f r0 = A_.row(0); + const double pr0 = dot(r0, p); + const double r0sqr = dot(r0, r0); + const double lhs = sqr(pr0); + const double rhs = r0sqr * distSqr * SQR_COS_ALPHA; + return lhs < rhs; + } + case 2: { + const Vec3f r0 = A_.row(0); + const Vec3f r1 = A_.row(1); + const Vec3f n = cross(r0, r1); + const double np = dot(n, p); + const double nsqr = dot(n, n); + const double lhs = sqr(np); + const double rhs = nsqr * distSqr * SQR_SIN_ALPHA; + return lhs > rhs; + } + default: + throw; + } + } +}; + +class MemorylessDecimator { + using Mesh = TriangleMesh; + +public: + float cost(const Mesh& mesh, const Graph::CollapseContext& context) const { + const Vec3f p0 = mesh.point(context.remaining); + const Vec3f p1 = mesh.point(context.removed); + const double lengthSqr = normSqr(p0 - p1); + std::set triangles = getTriangles(mesh, context); + std::set ring = getVertexRing(mesh, context); + std::set boundary; //= getBoundary(mesh, context); + Vec3f p = computePlacement(mesh, triangles, ring, boundary, lengthSqr); + /*if (!p || !isCollapseAllowed(ci, *p)) { + return Base::ILLEGAL_COLLAPSE; + }*/ + + const double volumeCost = computeVolumeCost(mesh, triangles, p); + const double boundaryCost = computeBoundaryCost(mesh, boundary, p); + + const double totalCost = 0.5 * boundaryCost * lengthSqr + 0.5 * volumeCost; + PVL_ASSERT(totalCost >= 0.); + + return totalCost; + } + + Vec3f placement(const Mesh& mesh, const Graph::CollapseContext& context) const { + const Vec3f p0 = mesh.point(context.remaining); + const Vec3f p1 = mesh.point(context.removed); + const double lengthSqr = normSqr(p0 - p1); + std::set faces = getTriangles(mesh, context); + std::set ring = getVertexRing(mesh, context); + std::set boundary; //= getBoundary(ci); + return computePlacement(mesh, faces, ring, boundary, lengthSqr); + } + + void postprocess(const Mesh&, const Graph::CollapseContext&) {} + +private: + inline Mat33f crossProductMatrixSqr(const Vec3f& p) const { + return { + Vec3f(p[1] * p[1] + p[2] * p[2], -p[0] * p[1], -p[0] * p[2]), + Vec3f(-p[0] * p[1], p[0] * p[0] + p[2] * p[2], -p[1] * p[2]), + Vec3f(-p[0] * p[2], -p[1] * p[2], p[0] * p[0] + p[1] * p[1]), + }; + } + + Vec3f computePlacement(const Mesh& mesh, + const std::set& triangles, + const std::set& ring, + const std::set& boundary, + const double lengthSqr) const { + LindstromTurkConstraints constraints; + + addVolumeAndBoundaryOptimizationConstraint( + mesh, constraints, triangles, boundary, lengthSqr); + std::cout << "after optimization - " << constraints.count() << " constraints" + << std::endl; + addBoundaryConstraint(mesh, constraints, boundary); + addVolumeConstraint(mesh, constraints, triangles); + std::cout << "after volume - " << constraints.count() << " constraints" << std::endl; + addShapeConstraint(mesh, constraints, ring); + std::cout << "after shape - " << constraints.count() << " constraints" << std::endl; + + Vec3f p = constraints.getPlacement(); + return p; + } + + void addVolumeConstraint(const Mesh& mesh, + LindstromTurkConstraints& constraints, + const std::set& faces) const { + if (constraints.full()) { + return; + } + Vec3f sumN(0); + double sumL = 0.; + + for (FaceHandle fh : faces) { + std::array tri = mesh.triangle(fh); + const Vec3f n = cross(tri[1] - tri[0], tri[2] - tri[0]); + const double l = dot(cross(tri[0], tri[1]), tri[2]); + sumN += n; + sumL += l; + } + + constraints.addConstraint(sumN, sumL); + } + + void addBoundaryConstraint(const Mesh& mesh, + LindstromTurkConstraints& constraints, + const std::set& edges) const { + if (edges.empty() || constraints.full()) { + return; + } + + Vec3f e1(0); + Vec3f e2(0); + + for (HalfEdgeHandle eh : edges) { + const Vec3f from = mesh.point(mesh.from(eh)); + const Vec3f to = mesh.point(mesh.to(eh)); + e1 += from - to; + e2 += cross(from, to); + } + + const Mat33f H = crossProductMatrixSqr(e1); + const Vec3f c = cross(e1, e2); + + constraints.addQuadraticConstraint(H, c); + } + + void addShapeConstraint(const Mesh& mesh, + LindstromTurkConstraints& constraints, + const std::set& ring) const { + if (constraints.full()) { + return; + } + Mat33f H = Mat33f::identity() * double(ring.size()); + Vec3f c(0); + for (VertexHandle v : ring) { + c -= mesh.point(v); + } + constraints.addQuadraticConstraint(H, c); + } + + void addVolumeAndBoundaryOptimizationConstraint(const Mesh& mesh, + LindstromTurkConstraints& constraints, + const std::set& faces, + const std::set& boundary, + const double lengthSqr) const { + if (constraints.full()) { + std::cout << "Skipping optimization contraint" << std::endl; + return; + } + + Mat33f H = Mat33f::null(); + Vec3f c(0, 0, 0); + + for (FaceHandle face : faces) { + std::array tri = mesh.triangle(face); + const Vec3f n = cross(tri[1] - tri[0], tri[2] - tri[0]); + const double l = dot(cross(tri[0], tri[1]), tri[2]); + + H += outerProd(n, n); + c -= n * l; + } + + if (!boundary.empty()) { + Mat33f Hb = Mat33f::null(); + Vec3f cb(0, 0, 0); + for (HalfEdgeHandle edge : boundary) { + const Vec3f from = mesh.point(mesh.from(edge)); + const Vec3f to = mesh.point(mesh.to(edge)); + const Vec3f n = cross(from, to); + const Vec3f dir = from - to; + + Hb += crossProductMatrixSqr(dir); + cb += cross(dir, n); + } + // 9 * boundary weight * homogenizing factor + const double volumeWeight = 0.5; + const double boundaryWeight = 0.5; + const double w = 9 * boundaryWeight * lengthSqr; + + H = H * volumeWeight + Hb * w; + c = c * volumeWeight + cb * w; + } + constraints.addQuadraticConstraint(H, c); + } + + std::set getTriangles(const Mesh& mesh, + const Graph::CollapseContext& context) const { + std::set faces; + for (FaceHandle fh : mesh.faceRing(context.remaining)) { + faces.insert(fh); + } + for (FaceHandle fh : mesh.faceRing(context.removed)) { + faces.insert(fh); + } + faces.erase(context.left); + faces.erase(context.right); + return faces; + } + + /*std::set getBoundary(const CollapseInfo& ci) const { + std::set handles; + for (VertexHandle v : { ci.v0, ci.v1 }) { + // incomming halfedges + for (auto iter = mesh_.vih_iter(v); iter.is_valid(); ++iter) { + HalfedgeHandle eh1 = *iter; + HalfedgeHandle eh2 = mesh_.opposite_halfedge_handle(eh1); + if (mesh_.is_boundary(eh1)) { + handles.insert(eh1); + } + if (mesh_.is_boundary(eh2)) { + handles.insert(eh2); + } + } + } + return handles; + }*/ + + std::set getVertexRing(const Mesh& mesh, + const Graph::CollapseContext& context) const { + std::set ring; + for (VertexHandle v1 : { context.removed, context.remaining }) { + // incomming halfedges + for (VertexHandle v2 : mesh.vertexRing(v1)) { + if (v2 != context.removed && v2 != context.remaining) { + ring.insert(v2); + } + } + } + return ring; + } + + double computeVolumeCost(const Mesh& mesh, + const std::set& faces, + const Vec3f& p) const { + double cost = 0.; + for (FaceHandle face : faces) { + std::array tri = mesh.triangle(face); + + const Vec3f v01 = tri[1] - tri[0]; + const Vec3f v02 = tri[2] - tri[0]; + const Vec3f n = cross(v01, v02); + const double l = dot(cross(tri[0], tri[1]), tri[2]); + cost += sqr(dot(n, p) - l); + } + return cost / 36.; + } + + double computeBoundaryCost(const Mesh& mesh, + const std::set& edges, + const Vec3f& p) const { + double cost = 0.; + for (HalfEdgeHandle edge : edges) { + Vec3f from = mesh.point(mesh.from(edge)); + Vec3f to = mesh.point(mesh.to(edge)); + + const Vec3f dir = from - to; + const Vec3f c = cross(dir, from - p); + cost += dot(c, c); + } + return cost / 4.; + } + + /* bool isCollapseAllowed(const CollapseInfo& ci, const Point& p) const { + // simulate collapse + mesh_.set_point(ci.v0, p); + mesh_.set_point(ci.v1, p); + OnScopeExit guard([&] { + // rollback the collapse + mesh_.set_point(ci.v0, ci.p0); + mesh_.set_point(ci.v1, ci.p1); + }); + + for (VertexHandle vh : { ci.v0, ci.v1 }) { + for (auto iter = mesh_.vf_iter(vh); iter.is_valid(); ++iter) { + FaceHandle face = *iter; + // these faces will be removed, so no need to test them + if (face != ci.fl && face != ci.fr) { + typename MeshT::Normal n1 = mesh_.normal(face); + typename MeshT::Normal n2 = mesh_.calc_face_normal(face); + const double cosPhi = dot(n1, n2); + if (cosPhi < MIN_COS_PHI) { + return false; + } + } + } + } + return true; + }*/ +}; + +} // namespace Pvl diff --git a/PlyReader.hpp b/PlyReader.hpp index 1891d98..fe21a46 100644 --- a/PlyReader.hpp +++ b/PlyReader.hpp @@ -20,7 +20,8 @@ class PlyReader { std::string line; // std::string countTag("element vertex "); while (std::getline(in_, line)) { - /* if (line.size() > countTag.size() && line.substr(0, countTag.size()) == countTag) { + /* if (line.size() > countTag.size() && line.substr(0, countTag.size()) == + countTag) { }*/ if (line == "end_header") { @@ -42,15 +43,15 @@ class PlyReader { std::size_t numVertices = 0; std::size_t numFaces = 0; while (std::getline(in_, line)) { - sscanf(line.c_str(), "element vertex %d", &numVertices); - sscanf(line.c_str(), "element face %d", &numFaces); + sscanf(line.c_str(), "element vertex %zu", &numVertices); + sscanf(line.c_str(), "element face %zu", &numFaces); if (line == "end_header") { break; } } TriangleMesh mesh; - std::cout << "Loading mesh with " << numVertices << " vertices and " << numFaces << " faces" - << std::endl; + std::cout << "Loading mesh with " << numVertices << " vertices and " << numFaces + << " faces" << std::endl; for (std::size_t i = 0; i < numVertices; ++i) { std::getline(in_, line); std::stringstream ss(line); diff --git a/PlyWriter.hpp b/PlyWriter.hpp index 9496bf0..04b5a8f 100644 --- a/PlyWriter.hpp +++ b/PlyWriter.hpp @@ -71,12 +71,12 @@ class PlyWriter { template PlyWriter& operator<<(const TriangleMesh& mesh) { - for (Index i = 0; i < mesh.numVertices(); ++i) { + for (Index i = 0; i < Index(mesh.numVertices()); ++i) { const Vec& p = mesh.points[i]; out_ << p[0] << " " << p[1] << " " << p[2] << " 0 0 0\n"; } std::size_t validCnt = 0; - for (Index i = 0; i < mesh.numFaces(); ++i) { + for (Index i = 0; i < Index(mesh.numFaces()); ++i) { if (!mesh.valid(FaceHandle(i))) { continue; } diff --git a/QuadricDecimator.hpp b/QuadricDecimator.hpp new file mode 100644 index 0000000..3015522 --- /dev/null +++ b/QuadricDecimator.hpp @@ -0,0 +1,79 @@ +#pragma once + +#include "Matrix.hpp" +#include "TriangleMesh.hpp" +#include +#include +#include + +namespace Pvl { + +/// \todo make callable with both Context and EdgeHandle? (distinguish using traits) +class QuadricDecimator { + using Mesh = TriangleMesh; + + std::map quadrics_; + +public: + QuadricDecimator(Mesh& mesh) { + for (VertexHandle vh : mesh.vertexRange()) { + quadrics_[vh] = Mat44f::null(); + } + + for (FaceHandle fh : mesh.faceRange()) { + Vec3f n = mesh.normal(fh); + Vec3f p0 = mesh.triangle(fh)[0]; + Vec4f plane; + plane[0] = n[0]; + plane[1] = n[1]; + plane[2] = n[2]; + plane[3] = -dot(p0, n); + + /// \todo optimize - compute together with normal + const float area = mesh.area(fh); + PVL_ASSERT(area > 0.f); + + Mat44f Q = outerProd(plane, plane) * area; + + std::array idxs = mesh.faceIndices(fh); + for (int i : idxs) { + quadrics_[VertexHandle(i)] += Q; + } + } + } + + float cost(const Mesh& mesh, const Graph::CollapseContext& context) const { + /*Vec3f v1 = mesh.point(context.remaining); + Vec3f v2 = mesh.point(context.removed);*/ + Mat44f Q = quadrics_.at(context.remaining) + quadrics_.at(context.removed); + Vec3f target = placement(mesh, context); + Vec4f v = homogeneous(target); + float result = dot(Q.transform(v), v); + // PVL_ASSERT(result >= 0.f); + return result; + } + + Vec3f placement(const Mesh& mesh, const Graph::CollapseContext& context) const { + Mat44f Q = quadrics_.at(context.remaining) + quadrics_.at(context.removed); + float eps = 1.e-6 * pow<4>(trace(Q)); + + Q(0, 3) = Q(1, 3) = Q(2, 3) = 0; + Q(3, 3) = 1; + // bit of magic to get dimensionless eps + float det = std::abs(determinant(Q)); + if (det > eps) { + Mat44f Qinv = invert(Q); + return Vec3f(Qinv(3, 0), Qinv(3, 1), Qinv(3, 2)); + } else { + Vec3f v1 = mesh.point(context.remaining); + Vec3f v2 = mesh.point(context.removed); + return 0.5 * (v1 + v2); + } + } + + void postprocess(const Mesh&, const Graph::CollapseContext& context) { + quadrics_[context.remaining] += quadrics_[context.removed]; + } +}; + +} // namespace Pvl diff --git a/Simplification.hpp b/Simplification.hpp index 1913ac3..8f062e8 100644 --- a/Simplification.hpp +++ b/Simplification.hpp @@ -7,19 +7,23 @@ namespace Pvl { template -class MidpointPlacement { -public: - Vec operator()(const Vec& v1, const Vec& v2) const { - return 0.5 * (v1 + v2); +struct SimpleDecimator { + template + typename Vec::Float cost(const Mesh& mesh, const Graph::CollapseContext& context) const { + Vec v1 = mesh.point(context.remaining); + Vec v2 = mesh.point(context.removed); + return normSqr(v1 - v2); } -}; -template -class EdgeLengthCost { -public: - typename Vec::Float operator()(const Vec& v1, const Vec& v2) const { - return norm(v1 - v2); + template + Vec placement(const Mesh& mesh, const Graph::CollapseContext& context) const { + Vec v1 = mesh.point(context.remaining); + Vec v2 = mesh.point(context.removed); + return 0.5 * (v1 + v2); } + + template + void postprocess(const Mesh&, const Graph::CollapseContext&) const {} }; class EdgeCountStop { @@ -35,43 +39,47 @@ class EdgeCountStop { }; class CollapseQueue { - using EdgeCost = std::pair; - std::map edges_; + using EdgeCost = std::pair; + std::map edges_; std::set queue_; public: - void insert(HalfEdgeHandle eh, const float cost) { + void insert(EdgeHandle eh, const float cost) { EdgeCost ec{ cost, eh }; PVL_ASSERT(edges_.find(eh) == edges_.end()); queue_.insert(ec); - edges_[eh] = ec; + edges_.insert(std::make_pair(eh, ec)); } - void update(HalfEdgeHandle eh, const float cost) { - /// \todo fix PVL_ASSERT(edges_.find(eh) != edges_.end()); - if (edges_.find(eh) == edges_.end()) { - insert(eh, cost); - } else { - EdgeCost oldEc = edges_[eh]; - queue_.erase(oldEc); + /*void update(EdgeHandle eh, const float cost) { + PVL_ASSERT(edges_.find(eh) != edges_.end()); + EdgeCost oldEc = edges_[eh]; + queue_.erase(oldEc); - EdgeCost newEc{ cost, eh }; - queue_.insert(newEc); - } - } + PVL_ASSERT(oldEc.second == eh); + EdgeCost newEc{ cost, eh }; + queue_.insert(newEc); + edges_[eh] = newEc; + }*/ - std::pair pop() { + std::pair pop() { EdgeCost ec = *queue_.begin(); queue_.erase(ec); - edges_.erase(ec.second); + std::size_t erased = edges_.erase(ec.second); + PVL_ASSERT(erased == 1); return std::make_pair(ec.second, ec.first); } - void remove(HalfEdgeHandle eh) { + void remove(EdgeHandle eh) { + /// \todo + if (edges_.find(eh) == edges_.end()) { + return; + } PVL_ASSERT(edges_.find(eh) != edges_.end()); EdgeCost ec = edges_[eh]; queue_.erase(ec); - edges_.erase(ec.second); + PVL_ASSERT(eh == ec.second); + edges_.erase(eh); } bool empty() const { @@ -139,46 +147,59 @@ void savePatch(const TriangleMesh& mesh, HalfEdgeHandle eh) { } } -template -void simplify(TriangleMesh& mesh, - const Cost& cost, - const Placement& placement, - const Stop& stop) { - +/// \todo add stop to decimator? (cost stop) +template +void simplify(TriangleMesh& mesh, Decimator& decimator, const Stop& stop) { CollapseQueue queue; for (EdgeHandle eh : mesh.edgeRange()) { - HalfEdgeHandle heh = mesh.halfEdge(eh); - float c = cost(mesh.points[mesh.from(heh)], mesh.points[mesh.to(heh)]); - queue.insert(heh, c); + float cost = decimator.cost(mesh, Graph::CollapseContext(mesh, eh)); + queue.insert(eh, cost); } /*for (HalfEdgeHandle heh : mesh.halfEdgeRange()) { + if (!mesh.valid(heh)) { + continue; + } EdgeHandle eh = mesh.edge(heh); HalfEdgeHandle h = mesh.halfEdge(eh); - float c = cost(mesh.points[mesh.from(h)], mesh.points[mesh.to(h)]); - queue.update(h, c); + float cost = decimator.cost(mesh, Graph::CollapseContext(mesh, heh)); + queue.update(h, cost); }*/ std::size_t cnt = 0; - HalfEdgeHandle collapsedEdge; + EdgeHandle collapsedEdge; float c; while (!queue.empty()) { std::tie(collapsedEdge, c) = queue.pop(); if (mesh.removed(collapsedEdge)) { continue; } - if (mesh.collapseAllowed(collapsedEdge)) { - Vec p = placement( - mesh.points[mesh.from(collapsedEdge)], mesh.points[mesh.to(collapsedEdge)]); - std::vector ring; + Graph::CollapseContext context(mesh, collapsedEdge); + if (mesh.collapseAllowed(context.edge)) { + Vec target = decimator.placement(mesh, context); + /*std::vector ring; for (HalfEdgeHandle eh : mesh.halfEdgeRing(mesh.from(collapsedEdge))) { ring.push_back(eh); } for (HalfEdgeHandle eh : mesh.halfEdgeRing(mesh.to(collapsedEdge))) { ring.push_back(eh); } - std::cout << "# " << cnt << " Collapsing " << mesh.to(collapsedEdge) << " into " - << mesh.from(collapsedEdge) << ", cost = " << c << std::endl; + */ + std::set ring; + for (VertexHandle vh : mesh.vertexRing(context.remaining)) { + if (vh != context.removed) { + ring.insert(vh); + queue.remove(mesh.edge(vh, context.remaining)); + } + } + for (VertexHandle vh : mesh.vertexRing(context.removed)) { + if (vh != context.remaining) { + ring.insert(vh); + queue.remove(mesh.edge(vh, context.removed)); + } + } + std::cout << "# " << cnt << " Collapsing " << mesh.to(context.edge) << " into " + << mesh.from(context.edge) << ", cost = " << c << std::endl; /*std::cout << "ring:\n"; for (HalfEdgeHandle eh : ring) { std::cout << mesh.from(eh) << "-" << mesh.to(eh) << "\n"; @@ -190,18 +211,18 @@ void simplify(TriangleMesh& mesh, } std::cout << std::endl; savePatch(mesh, collapsedEdge);*/ - mesh.collapse(collapsedEdge, p); - for (HalfEdgeHandle heh : ring) { - /// \todo fix this - if (!mesh.valid(heh)) { - continue; - } - HalfEdgeHandle heh1 = mesh.halfEdge(mesh.edge(heh)); - if (!mesh.removed(heh1)) { - float c = cost(mesh.points[mesh.from(heh1)], mesh.points[mesh.to(heh1)]); - queue.update(heh1, c); - } + mesh.collapse(context.edge, target); + for (VertexHandle vh : ring) { + EdgeHandle eh = mesh.edge(context.remaining, vh); + PVL_ASSERT(!mesh.removed(eh)); + PVL_ASSERT(mesh.valid(eh)); + + // HalEdgeHandle heh1 = mesh.halfEdge(mesh.edge(heh)); + // PVL_ASSERT(!mesh.removed(heh1)); + float cost = decimator.cost(mesh, Graph::CollapseContext(mesh, eh)); + queue.insert(eh, cost); } + decimator.postprocess(mesh, context); if (stop(cnt++)) { return; } diff --git a/TriangleMesh.hpp b/TriangleMesh.hpp index 6802f5b..63cba79 100644 --- a/TriangleMesh.hpp +++ b/TriangleMesh.hpp @@ -17,6 +17,28 @@ class TriangleMesh : public Graph { std::vector points; + Vec3f point(VertexHandle vh) const { + return points[vh]; + } + + std::array triangle(FaceHandle fh) const { + std::array idxs = faceIndices(fh); + Vec3f p0 = points[idxs[0]]; + Vec3f p1 = points[idxs[1]]; + Vec3f p2 = points[idxs[2]]; + return { p0, p1, p2 }; + } + + Vec3f normal(FaceHandle fh) const { + std::array tr = triangle(fh); + return normalize(cross(tr[1] - tr[0], tr[2] - tr[0])); + } + + float area(FaceHandle fh) const { + std::array tr = triangle(fh); + return 0.5f * norm(cross(tr[1] - tr[0], tr[2] - tr[0])); + } + // to -> from void collapse(HalfEdgeHandle eh, const Vec& placement) { points[to(eh)] = points[from(eh)] = placement;