Skip to content

Commit

Permalink
added quadric and memoryless simplification
Browse files Browse the repository at this point in the history
  • Loading branch information
@sevec
sevec committed Jun 24, 2020
1 parent e3495cb commit acc3d87
Show file tree
Hide file tree
Showing 10 changed files with 770 additions and 129 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
32 changes: 19 additions & 13 deletions Graph.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
};

Expand All @@ -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) {
Expand Down Expand Up @@ -762,12 +766,14 @@ class Graph {
std::vector<VertexHandle> 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) {
Expand Down Expand Up @@ -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() {}
Expand Down
7 changes: 6 additions & 1 deletion Math.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,18 @@ struct Pow<T, 2> {
return value * value;
}
};

template <typename T>
struct Pow<T, 3> {
T operator()(const T& value) {
return value * value * value;
}
};
template <typename T>
struct Pow<T, 4> {
T operator()(const T& value) {
return sqr(value * value);
}
};

template <int N, typename T>
T pow(const T& value) {
Expand Down
175 changes: 125 additions & 50 deletions Matrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,12 +32,16 @@ class Matrix {
return rows_[r][c];
}


const Vector<T, Cols>& row(const int idx) const {
PVL_ASSERT(idx < Rows);
return rows_[idx];
}

Vector<T, Cols>& row(const int idx) {
PVL_ASSERT(idx < Rows);
return rows_[idx];
}

Vector<T, Rows> column(const int idx) const {
PVL_ASSERT(idx < Cols);
Vector<T, Rows> col;
Expand Down Expand Up @@ -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;
Expand All @@ -114,16 +139,26 @@ class Matrix {
using Mat33f = Matrix<float, 3, 3>;
using Mat44f = Matrix<float, 4, 4>;

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 <typename T, int N, int M>
Matrix<T, N, M> outerProd(const Vector<T, N>& v1, const Vector<T, M>& v2) {
Matrix<T, N, M> 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 <typename T, int N>
inline T trace(const Matrix<T, N, N>& 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)) +
Expand All @@ -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++) {
Expand Down
Loading

0 comments on commit acc3d87

Please sign in to comment.