From 127033e6ba1d66164e2387a1c3128a89a375967d Mon Sep 17 00:00:00 2001
From: pavel <sevecek@sirrah.troja.mff.cuni.cz>
Date: Tue, 30 Jun 2020 14:11:11 +0200
Subject: [PATCH] changed default index type to uint32_t

---
 Assert.hpp           |   4 ++
 Graph.hpp            |  30 +++-----
 Optional.hpp         |   1 +
 PlyReader.hpp        |   8 ++-
 PlyWriter.hpp        |   8 +--
 QuadricDecimator.hpp |  10 +--
 Refinement.hpp       |  12 ++--
 Simplification.hpp   |  31 ++++++--
 TriangleMesh.hpp     |  15 ++--
 Vector.hpp           |   5 +-
 main.cpp             | 165 ++++++++++++++++++++++++-------------------
 11 files changed, 170 insertions(+), 119 deletions(-)

diff --git a/Assert.hpp b/Assert.hpp
index 75b2df8..7058a7e 100644
--- a/Assert.hpp
+++ b/Assert.hpp
@@ -2,4 +2,8 @@
 
 #include <cassert>
 
+#ifndef NDEBUG
 #define PVL_ASSERT(x, ...) assert(x)
+#else
+#define PVL_ASSERT(x, ...) (void)sizeof(x)
+#endif
diff --git a/Graph.hpp b/Graph.hpp
index 5fdc3b8..e1065a0 100644
--- a/Graph.hpp
+++ b/Graph.hpp
@@ -14,8 +14,10 @@
 
 namespace Pvl {
 
-template <typename TObject, typename Index = int>
+template <typename TObject>
 class Handle {
+    using Index = std::uint32_t;
+
     Index idx_;
 
 public:
@@ -67,7 +69,7 @@ struct HalfEdge {
     }
 
     bool boundary() const {
-        return opposite < 0;
+        return opposite == HalfEdgeHandle(-1);
     }
 };
 
@@ -288,10 +290,6 @@ class Graph {
         return boundary(vertices_[vh].edge);
     }
     VertexHandle to(HalfEdgeHandle eh) const {
-        if (!valid(eh)) {
-            std::cout << "Invalid he handle " << halfEdges_[halfEdges_[eh].prev].to << "-"
-                      << halfEdges_[eh].to << std::endl;
-        }
         PVL_ASSERT(valid(eh));
         return halfEdges_[eh].to;
     }
@@ -379,19 +377,19 @@ class Graph {
         return prev(emanating(vh));
     }
     bool valid(FaceHandle fh) const {
-        PVL_ASSERT(fh < int(faces_.size()));
+        PVL_ASSERT(fh < faces_.size());
         return fh != FaceHandle(-1) && faces_[fh].edge != HalfEdgeHandle(-1);
     }
     bool valid(HalfEdgeHandle heh) const {
-        PVL_ASSERT(heh < int(halfEdges_.size()));
+        PVL_ASSERT(heh < halfEdges_.size());
         return heh != HalfEdgeHandle(-1) && halfEdges_[heh].left != FaceHandle(-1);
     }
     bool valid(VertexHandle vh) const {
-        PVL_ASSERT(vh < int(vertices_.size()));
+        PVL_ASSERT(vh < vertices_.size());
         return vh != VertexHandle(-1) && vertices_[vh].edge != HalfEdgeHandle(-1);
     }
     bool valid(EdgeHandle eh) const {
-        PVL_ASSERT(eh < int(halfEdges_.size()));
+        PVL_ASSERT(eh < halfEdges_.size());
         HalfEdgeHandle heh(eh.index());
         return valid(heh) && (boundary(heh) || from(heh) < to(heh));
     }
@@ -469,7 +467,7 @@ class Graph {
 
     private:
         bool end() const {
-            return eh_ >= int(graph_.halfEdges_.size());
+            return eh_ >= graph_.halfEdges_.size();
         }
 
         bool dereferencable() const {
@@ -712,28 +710,22 @@ class Graph {
             PVL_ASSERT(!onesided);
             ev0 = opposite(prev(heh));
         }
-        std::cout << "Assigning emanating halfedge for v0 - " << from(ev0) << "-" << to(ev0)
-                  << std::endl;
         PVL_ASSERT(valid(ev0));
 
         if (vertices_[vL].edge == prev(heh)) {
             // move to any other halfedge, it cannot be a boundary
             PVL_ASSERT(!boundary(prev(heh)));
             evL = next(opposite(prev(heh)));
-            std::cout << "Assigning emanating halfedge for evL - " << from(evL) << "-"
-                      << to(evL) << std::endl;
         } else {
             // keep the edge as it might be boundary
             evL = vertices_[vL].edge;
         }
-        PVL_ASSERT(evL);
+        PVL_ASSERT(valid(evL));
 
         if (!onesided) {
             if (vertices_[vR].edge == prev(oheh)) {
                 PVL_ASSERT(!boundary(prev(oheh)));
                 evR = opposite(next(oheh));
-                std::cout << "Assigning emanating halfedge for evR - " << from(evR) << "-"
-                          << to(evR) << std::endl;
             } else {
                 evR = vertices_[vR].edge;
             }
@@ -881,7 +873,7 @@ class Graph {
     }*/
 
     bool removed(HalfEdgeHandle heh) {
-        return halfEdges_[heh].left == -1;
+        return halfEdges_[heh].left == HalfEdgeHandle(-1);
     }
 
     bool removed(EdgeHandle eh) {
diff --git a/Optional.hpp b/Optional.hpp
index b37dbf4..97d1dbf 100644
--- a/Optional.hpp
+++ b/Optional.hpp
@@ -1,6 +1,7 @@
 #pragma once
 #include "Assert.hpp"
 #include <type_traits>
+#include <utility>
 
 namespace Pvl {
 
diff --git a/PlyReader.hpp b/PlyReader.hpp
index fe21a46..3d91175 100644
--- a/PlyReader.hpp
+++ b/PlyReader.hpp
@@ -38,7 +38,7 @@ class PlyReader {
         return points;
     }
 
-    TriangleMesh<Vec3f, int> readMesh() {
+    TriangleMesh<Vec3f> readMesh() {
         std::string line;
         std::size_t numVertices = 0;
         std::size_t numFaces = 0;
@@ -49,7 +49,7 @@ class PlyReader {
                 break;
             }
         }
-        TriangleMesh<Vec3f, int> mesh;
+        TriangleMesh<Vec3f> mesh;
         std::cout << "Loading mesh with " << numVertices << " vertices and " << numFaces
                   << " faces" << std::endl;
         for (std::size_t i = 0; i < numVertices; ++i) {
@@ -60,12 +60,16 @@ class PlyReader {
             mesh.addVertex();
             mesh.points.push_back(p);
         }
+        std::cout << "Added " << mesh.numVertices() << " vertices " << std::endl;
         for (std::size_t i = 0; i < numFaces; ++i) {
             std::getline(in_, line);
             std::stringstream ss(line);
             int dummy, a, b, c;
             ss >> dummy >> a >> b >> c;
             mesh.addFace(VertexHandle(a), VertexHandle(b), VertexHandle(c));
+            if (i % 100000 == 99999) {
+                std::cout << "added " << mesh.numFaces() << " faces" << std::endl;
+            }
         }
         return mesh;
     }
diff --git a/PlyWriter.hpp b/PlyWriter.hpp
index a6f2fde..f61ad16 100644
--- a/PlyWriter.hpp
+++ b/PlyWriter.hpp
@@ -69,14 +69,14 @@ class PlyWriter {
         return *this;
     }
 
-    template <typename Vec, typename Index>
-    PlyWriter& operator<<(const TriangleMesh<Vec, Index>& mesh) {
-        for (Index i = 0; i < Index(mesh.numVertices()); ++i) {
+    template <typename Vec>
+    PlyWriter& operator<<(const TriangleMesh<Vec>& mesh) {
+        for (std::uint32_t i = 0; i < 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 < Index(mesh.numFaces()); ++i) {
+        for (std::uint32_t i = 0; i < mesh.numFaces(); ++i) {
             if (!mesh.valid(FaceHandle(i))) {
                 continue;
             }
diff --git a/QuadricDecimator.hpp b/QuadricDecimator.hpp
index 737b47a..1ff2a00 100644
--- a/QuadricDecimator.hpp
+++ b/QuadricDecimator.hpp
@@ -28,7 +28,12 @@ class QuadricDecimator {
         }
 
         for (FaceHandle fh : mesh.faceRange()) {
-            Point n = mesh.normal(fh);
+            Point n = mesh.areaNormal(fh);
+            const Float area = norm(n);
+            if (area == 0) {
+                continue;
+            }
+            n /= area;
             Point p0 = mesh.triangle(fh)[0]; /// \todo fix
             Plane plane;
             plane[0] = n[0];
@@ -36,9 +41,6 @@ class QuadricDecimator {
             plane[2] = n[2];
             plane[3] = -dotProd(p0, n);
 
-            /// \todo optimize - compute together with normal
-            const Float area = mesh.area(fh);
-            PVL_ASSERT(area > 0);
 
             Quadric Q = outerProd(plane, plane) * area;
 
diff --git a/Refinement.hpp b/Refinement.hpp
index c014407..9284d14 100644
--- a/Refinement.hpp
+++ b/Refinement.hpp
@@ -6,8 +6,10 @@
 
 namespace Pvl {
 
-template <typename ConcurrencyTag = SequentialTag, typename Vec, typename Index>
-void laplacianSmoothing(TriangleMesh<Vec, Index>& mesh, bool preserveBoundary = true, float rer = 1.f) {
+template <typename ConcurrencyTag = SequentialTag, typename Vec>
+void laplacianSmoothing(TriangleMesh<Vec>& mesh,
+    bool preserveBoundary = true,
+    float rer = 1.f) {
     std::vector<Vec> laplacian(mesh.numVertices(), Vec(0));
     ParallelForEach<ConcurrencyTag>()(
         mesh.vertexRange(), [&mesh, &laplacian, preserveBoundary](VertexHandle v1) {
@@ -26,8 +28,8 @@ void laplacianSmoothing(TriangleMesh<Vec, Index>& mesh, bool preserveBoundary =
         });
     std::vector<Vec> biharmonic(mesh.numVertices(), Vec(0));
     if (rer > 0.f) {
-        ParallelForEach<ConcurrencyTag>()(
-            mesh.vertexRange(), [&mesh, &laplacian, &biharmonic, preserveBoundary](VertexHandle v1) {
+        ParallelForEach<ConcurrencyTag>()(mesh.vertexRange(),
+            [&mesh, &laplacian, &biharmonic, preserveBoundary](VertexHandle v1) {
                 if (preserveBoundary && mesh.boundary(v1)) {
                     return;
                 }
@@ -46,7 +48,7 @@ void laplacianSmoothing(TriangleMesh<Vec, Index>& mesh, bool preserveBoundary =
     }
 
     ParallelFor<ConcurrencyTag>()(
-        Index(0), Index(mesh.numVertices()), [&mesh, &laplacian, &biharmonic, &rer](std::size_t i) {
+        0, mesh.numVertices(), [&mesh, &laplacian, &biharmonic, &rer](std::size_t i) {
             mesh.points[i] += 0.5 * (rer * biharmonic[i] + (1.f - rer) * laplacian[i]);
         });
 }
diff --git a/Simplification.hpp b/Simplification.hpp
index d8468bc..1cf9466 100644
--- a/Simplification.hpp
+++ b/Simplification.hpp
@@ -80,11 +80,24 @@ class EdgeCountStop {
     EdgeCountStop(std::size_t count)
         : count_(count) {}
 
-    bool operator()(std::size_t collapsed) const {
-        return collapsed > count_;
+    bool operator()(std::size_t, std::size_t collapsedEdges) const {
+        return collapsedEdges > count_;
     }
 };
 
+class FaceCountStop {
+    std::size_t count_;
+
+public:
+    FaceCountStop(std::size_t count)
+        : count_(count) {}
+
+    bool operator()(std::size_t collapsedFaces, std::size_t) const {
+        return collapsedFaces > count_;
+    }
+};
+
+
 class CollapseQueue {
     using EdgeCost = std::pair<float, EdgeHandle>;
     std::map<EdgeHandle, EdgeCost> edges_;
@@ -196,8 +209,8 @@ void savePatch(const TriangleMesh<Vec, Index>& mesh, HalfEdgeHandle eh) {
 }*/
 
 /// \todo add stop to decimator? (cost stop)
-template <typename Vec, typename Index, typename Decimator, typename Stop>
-void simplify(TriangleMesh<Vec, Index>& mesh, Decimator& decimator, const Stop& stop) {
+template <typename Vec, typename Decimator, typename Stop>
+void simplify(TriangleMesh<Vec>& mesh, Decimator& decimator, const Stop& stop) {
     using Float = typename Vec::Float;
 
     CollapseQueue queue;
@@ -208,7 +221,8 @@ void simplify(TriangleMesh<Vec, Index>& mesh, Decimator& decimator, const Stop&
         }
     }
 
-    std::size_t cnt = 0;
+    std::size_t numCollapsedEdges = 0;
+    std::size_t numCollapsedFaces = 0;
     EdgeHandle collapsedEdge;
     float c;
     while (!queue.empty()) {
@@ -271,7 +285,12 @@ void simplify(TriangleMesh<Vec, Index>& mesh, Decimator& decimator, const Stop&
                 }
             }
             decimator.postprocess(mesh, context);
-            if (stop(cnt++)) {
+            if (numCollapsedEdges % 10000 == 0) {
+                std::cout << "# " << numCollapsedEdges << " collapsed" << std::endl;
+            }
+            numCollapsedEdges++;
+            numCollapsedFaces += 1 + int(context.right != FaceHandle(-1));
+            if (stop(numCollapsedFaces, numCollapsedEdges)) {
                 return;
             }
         }
diff --git a/TriangleMesh.hpp b/TriangleMesh.hpp
index 4c4ae3c..2788a68 100644
--- a/TriangleMesh.hpp
+++ b/TriangleMesh.hpp
@@ -6,7 +6,7 @@
 
 namespace Pvl {
 
-template <typename Vec, typename Index>
+template <typename Vec>
 class TriangleMesh : public Graph {
 public:
     using Point = Vec;
@@ -30,14 +30,19 @@ class TriangleMesh : public Graph {
         return { p0, p1, p2 };
     }
 
-    Point normal(FaceHandle fh) const {
+    // normal scaled by area
+    Point areaNormal(FaceHandle fh) const {
         std::array<Point, 3> tr = triangle(fh);
-        return normalize(crossProd(tr[1] - tr[0], tr[2] - tr[0]));
+        return 0.5 * crossProd(tr[1] - tr[0], tr[2] - tr[0]);
+    }
+
+    // normalized normal
+    Point normal(FaceHandle fh) const {
+        return normalize(areaNormal(fh));
     }
 
     float area(FaceHandle fh) const {
-        std::array<Vec3f, 3> tr = triangle(fh);
-        return 0.5f * norm(crossProd(tr[1] - tr[0], tr[2] - tr[0]));
+        return norm(areaNormal(fh));
     }
 
     // to -> from
diff --git a/Vector.hpp b/Vector.hpp
index dfcbae1..34d2d2c 100644
--- a/Vector.hpp
+++ b/Vector.hpp
@@ -106,6 +106,7 @@ class Vector {
     }
 
     Vector operator/(const T f) const {
+        PVL_ASSERT(f != 0);
         Vector res;
         for (int i = 0; i < Dim; ++i) {
             res[i] = values_[i] / f;
@@ -180,7 +181,9 @@ T normL1(const Vector<T, Dim>& v) {
 
 template <typename T, int Dim>
 Vector<T, Dim> normalize(const Vector<T, Dim>& v) {
-    return v / norm(v);
+    T length = norm(v);
+    PVL_ASSERT(length > 0);
+    return v / length;
 }
 
 inline Vec3f crossProd(const Vec3f& v1, const Vec3f& v2) {
diff --git a/main.cpp b/main.cpp
index ff1eb97..a1be7d0 100644
--- a/main.cpp
+++ b/main.cpp
@@ -19,8 +19,8 @@
 
 using namespace Pvl;
 /*std::cout << "Test" << std::endl;
-Pvl::KdTree<Pvl::Vec3f> tree;
-std::vector<Pvl::Vec3f> points = {
+KdTree<Vec3f> tree;
+std::vector<Vec3f> points = {
     { 0.f, 0.f, 0.f },
     { 0.f, 0.f, 4.f },
     { 1.f, 0.f, 0.f },
@@ -29,7 +29,7 @@ tree.build(points);
 
 // std::vector<uint32_t> neighs;
 std::array<uint32_t, 3> neighs;
-int n = tree.rangeQuery(Pvl::Vec3f(0.f, 0.f, 0.f), 1.5f, neighs.begin());
+int n = tree.rangeQuery(Vec3f(0.f, 0.f, 0.f), 1.5f, neighs.begin());
 
 std::cout << "Found " << n << " neighs" << std::endl;
 for (auto& n : neighs) {
@@ -37,7 +37,7 @@ for (auto& n : neighs) {
 }
 
 std::stringstream ss;
-Pvl::PlyWriter ply(ss);
+PlyWriter ply(ss);
 for (auto& p : points) {
     ply << p;
 }
@@ -45,37 +45,37 @@ ply.close();
 std::cout << ss.str() << "\n";
 
 std::stringstream iss(ss.str());
-Pvl::PlyReader rd(iss);
-std::vector<Pvl::Vec3f> ps = rd.read();
+PlyReader rd(iss);
+std::vector<Vec3f> ps = rd.read();
 std::cout << "Read points \n";
 for (auto& p : ps) {
     std::cout << p[0] << "," << p[1] << "," << p[2] << "\n";
 }
 
 std::cout << "Grid\n";
-Pvl::UniformGrid<float, 3> grid(Pvl::Vec3i(2, 3, 2));
+UniformGrid<float, 3> grid(Vec3i(2, 3, 2));
 int cntr = 0;
 for (auto& f : grid) {
     std::cout << cntr++ << "-" << f << "\n";
 }*/
 
-/* Pvl::Graph graph;
- Pvl::VertexHandle a = graph.addVertex();
- Pvl::VertexHandle b = graph.addVertex();
- Pvl::VertexHandle c = graph.addVertex();
- Pvl::VertexHandle d = graph.addVertex();
- Pvl::VertexHandle e = graph.addVertex();
- Pvl::VertexHandle f = graph.addVertex();
+/* Graph graph;
+ VertexHandle a = graph.addVertex();
+ VertexHandle b = graph.addVertex();
+ VertexHandle c = graph.addVertex();
+ VertexHandle d = graph.addVertex();
+ VertexHandle e = graph.addVertex();
+ VertexHandle f = graph.addVertex();
 
- Pvl::FaceHandle f1 = graph.addFace(a, b, c);
- Pvl::FaceHandle f2 = graph.addFace(e, f, a);
- Pvl::FaceHandle f3 = graph.addFace(d, e, a);
- Pvl::FaceHandle f4 = graph.addFace(a, c, d);
+ FaceHandle f1 = graph.addFace(a, b, c);
+ FaceHandle f2 = graph.addFace(e, f, a);
+ FaceHandle f3 = graph.addFace(d, e, a);
+ FaceHandle f4 = graph.addFace(a, c, d);
 
  std::cout << "Faces = " << graph.numFaces() << std::endl;
  std::cout << "Vertices = " << graph.numVertices() << std::endl;
 
- for (Pvl::VertexHandle vh : graph.vertexRing(a)) {
+ for (VertexHandle vh : graph.vertexRing(a)) {
      std::cout << " - " << vh.index() << std::endl;
  }*/
 
@@ -208,7 +208,7 @@ TEST_CASE_METHOD(SimpleGraphFixture, "edge range", "[graph]") {
 }
 
 TEST_CASE_METHOD(SimpleGraphFixture, "collapse inner edge", "[simplify]") {
-    Graph::EdgeRange edges = graph.edgeRange();
+    // Graph::EdgeRange edges = graph.edgeRange();
     EdgeHandle collapsible = graph.edge(vA, vB);
 
     /*REQUIRE(std::all_of(edges.begin(), edges.end(), [&](EdgeHandle eh) {
@@ -276,12 +276,12 @@ TEST_CASE_METHOD(SimpleGraphFixture, "collapse boundary edge", "[simplify]") {
 
 TEST_CASE("collapse simple 2", "[simplify]") {
     Graph graph;
-    std::array<Pvl::VertexHandle, 8> vs;
+    std::array<VertexHandle, 8> vs;
     for (int i = 0; i < 8; ++i) {
         vs[i] = graph.addVertex();
     }
-    Pvl::VertexHandle vA = graph.addVertex();
-    Pvl::VertexHandle vB = graph.addVertex();
+    VertexHandle vA = graph.addVertex();
+    VertexHandle vB = graph.addVertex();
 
     graph.addFace(vs[0], vA, vs[1]);
     graph.addFace(vA, vs[2], vs[1]);
@@ -334,24 +334,24 @@ TEST_CASE("collapse simple 2", "[simplify]") {
 }
 
 TEST_CASE("simplify simple", "[simplify]") {
-    Pvl::TriangleMesh<Pvl::Vec3f, std::size_t> mesh;
-    std::array<Pvl::VertexHandle, 8> vs;
+    TriangleMesh<Vec3f> mesh;
+    std::array<VertexHandle, 8> vs;
     for (int i = 0; i < 8; ++i) {
         vs[i] = mesh.addVertex();
     }
-    Pvl::VertexHandle vA = mesh.addVertex();
-    Pvl::VertexHandle vB = mesh.addVertex();
+    VertexHandle vA = mesh.addVertex();
+    VertexHandle vB = mesh.addVertex();
     mesh.points.resize(mesh.numVertices());
-    mesh.points[vA] = Pvl::Vec3f(0, 0, 0);
-    mesh.points[vB] = Pvl::Vec3f(1, 0, 0);
-    mesh.points[vs[0]] = Pvl::Vec3f(-0.5, 0, 0);
-    mesh.points[vs[1]] = Pvl::Vec3f(-0.1, 0.5, 0);
-    mesh.points[vs[2]] = Pvl::Vec3f(0.5, 0.6, 0);
-    mesh.points[vs[3]] = Pvl::Vec3f(1.3, 0.3, 0);
-    mesh.points[vs[4]] = Pvl::Vec3f(1.6, 0, 0);
-    mesh.points[vs[5]] = Pvl::Vec3f(1.3, -0.3, 0);
-    mesh.points[vs[6]] = Pvl::Vec3f(1, -0.8, 0);
-    mesh.points[vs[7]] = Pvl::Vec3f(0, -0.9, 0);
+    mesh.points[vA] = Vec3f(0, 0, 0);
+    mesh.points[vB] = Vec3f(1, 0, 0);
+    mesh.points[vs[0]] = Vec3f(-0.5, 0, 0);
+    mesh.points[vs[1]] = Vec3f(-0.1, 0.5, 0);
+    mesh.points[vs[2]] = Vec3f(0.5, 0.6, 0);
+    mesh.points[vs[3]] = Vec3f(1.3, 0.3, 0);
+    mesh.points[vs[4]] = Vec3f(1.6, 0, 0);
+    mesh.points[vs[5]] = Vec3f(1.3, -0.3, 0);
+    mesh.points[vs[6]] = Vec3f(1, -0.8, 0);
+    mesh.points[vs[7]] = Vec3f(0, -0.9, 0);
 
     mesh.addFace(vs[0], vA, vs[1]);
     mesh.addFace(vA, vs[2], vs[1]);
@@ -364,42 +364,42 @@ TEST_CASE("simplify simple", "[simplify]") {
     mesh.addFace(vs[6], vB, vA);
     mesh.addFace(vs[0], vs[7], vA);
 
-    Pvl::EdgeHandle eh = mesh.edge(vA, vB);
+    EdgeHandle eh = mesh.edge(vA, vB);
     std::cout << "A-B edge = " << eh << std::endl;
     {
         std::ofstream ofs("base.ply");
-        Pvl::PlyWriter writer(ofs);
+        PlyWriter writer(ofs);
         writer << mesh;
     }
     std::cout << "A ring:" << std::endl;
-    for (Pvl::VertexHandle vh : mesh.vertexRing(vA)) {
+    for (VertexHandle vh : mesh.vertexRing(vA)) {
         std::cout << vh << ",";
     }
     std::cout << std::endl;
     std::cout << "B ring:" << std::endl;
-    for (Pvl::VertexHandle vh : mesh.vertexRing(vB)) {
+    for (VertexHandle vh : mesh.vertexRing(vB)) {
         std::cout << vh << ",";
     }
     std::cout << std::endl;
 
-    mesh.collapse(eh, Pvl::Vec3f(0.5, 0., 0.));
+    mesh.collapse(eh, Vec3f(0.5, 0., 0.));
     std::cout << "After collapse A ring:" << std::endl;
     std::cout << "A ring:" << std::endl;
-    for (Pvl::VertexHandle vh : mesh.vertexRing(vA)) {
+    for (VertexHandle vh : mesh.vertexRing(vA)) {
         std::cout << vh << ",";
     }
     std::cout << std::endl;
 
     {
         std::ofstream ofs("simplified.ply");
-        Pvl::PlyWriter writer(ofs);
+        PlyWriter writer(ofs);
         writer << mesh;
     }
 }
 
 TEST_CASE("Test bunny", "[mesh]") {
     std::ifstream ifs("/home/pavel/projects/pvl/data/bunny-fixed.ply");
-    Pvl::PlyReader reader(ifs);
+    PlyReader reader(ifs);
     auto bunny = reader.readMesh();
 
     auto vertices = bunny.vertexRange();
@@ -416,19 +416,19 @@ TEST_CASE("Test bunny", "[mesh]") {
 
 TEST_CASE("simplify bunny", "[simplify]") {
     std::ifstream ifs("/home/pavel/projects/pvl/data/bunny-fixed.ply");
-    Pvl::PlyReader reader(ifs);
+    PlyReader reader(ifs);
     auto bunny = reader.readMesh();
 
     std::cout << "Simplifying bunny " << std::endl;
-    Pvl::SimpleDecimator<decltype(bunny)> decimator;
-    Pvl::simplify(bunny, decimator, Pvl::EdgeCountStop{ 32500 });
+    SimpleDecimator<decltype(bunny)> decimator;
+    simplify(bunny, decimator, EdgeCountStop{ 32500 });
     /*for (int cnt : { 6000, 6000, 6000, 6000, 6000, 3000 }) {
-        Pvl::SimpleDecimator<Pvl::Vec3f> decimator;
-        Pvl::simplify(bunny, decimator, Pvl::EdgeCountStop{ cnt });
+        SimpleDecimator<Vec3f> decimator;
+        simplify(bunny, decimator, EdgeCountStop{ cnt });
     }*/
     {
         std::ofstream ofs("simplified-simple.ply");
-        Pvl::PlyWriter writer(ofs);
+        PlyWriter writer(ofs);
         writer << bunny;
     }
 }
@@ -440,8 +440,8 @@ TEST_CASE("simplify bunny 2", "[simplify]") {
 
     std::cout << "Simplifying bunny " << std::endl;
 
-    PreventFaceFoldDecorator<QuadricDecimator<TriangleMesh<Vec3f, int>>> decimator(bunny);
-    simplify(bunny, decimator, Pvl::EdgeCountStop{ 32500 });
+    PreventFaceFoldDecorator<QuadricDecimator<TriangleMesh<Vec3f>>> decimator(bunny);
+    simplify(bunny, decimator, EdgeCountStop{ 32500 });
 
 
     {
@@ -453,27 +453,46 @@ TEST_CASE("simplify bunny 2", "[simplify]") {
 
 TEST_CASE("simplify bunny 3", "[simplify]") {
     std::ifstream ifs("/home/pavel/projects/pvl/data/bunny-fixed.ply");
-    Pvl::PlyReader reader(ifs);
+    PlyReader reader(ifs);
     auto bunny = reader.readMesh();
 
     std::cout << "Simplifying bunny " << std::endl;
 
-    Pvl::MemorylessDecimator<decltype(bunny)> decimator; //(bunny);
-    Pvl::simplify(bunny, decimator, Pvl::EdgeCountStop{ 32500 });
+    MemorylessDecimator<decltype(bunny)> decimator; //(bunny);
+    simplify(bunny, decimator, EdgeCountStop{ 32500 });
     /*for (int cnt : { 6000, 6000, 6000, 6000, 6000, 3000 }) {
-        Pvl::QuadricDecimator decimator(bunny);
-        Pvl::simplify(bunny, decimator, Pvl::EdgeCountStop{ cnt });
+        QuadricDecimator decimator(bunny);
+        simplify(bunny, decimator, EdgeCountStop{ cnt });
     }*/
 
     {
         std::ofstream ofs("simplified-lt.ply");
-        Pvl::PlyWriter writer(ofs);
+        PlyWriter writer(ofs);
+        writer << bunny;
+    }
+}
+
+
+TEST_CASE("simplify denver", "[simplify]") {
+    std::ifstream ifs("/home/pavel/projects/pvl/data/denver.ply");
+    PlyReader reader(ifs);
+    auto bunny = reader.readMesh();
+
+    std::cout << "Simplifying denver " << std::endl;
+
+    PreventFaceFoldDecorator<QuadricDecimator<TriangleMesh<Vec3f>>> decimator(bunny);
+    simplify(bunny, decimator, FaceCountStop{ bunny.numFaces() - 100000 });
+    std::cout << "Simplification complete" << std::endl;
+
+    {
+        std::ofstream ofs("simplified-denver.ply");
+        PlyWriter writer(ofs);
         writer << bunny;
     }
 }
 /*   std::ifstream ifs("/home/pavel/projects/pvl/data/pc.ply");
-   Pvl::PlyReader reader(ifs);
-   std::vector<Pvl::Vec3f> points;
+   PlyReader reader(ifs);
+   std::vector<Vec3f> points;
    {
        auto all = reader.readCloud();
        for (std::size_t i = 0; i < all.size(); i += 10) {
@@ -482,54 +501,54 @@ TEST_CASE("simplify bunny 3", "[simplify]") {
    }
    std::cout << "Loaded cloud with " << points.size() << " points " << std::endl;
    std::cout << "Estimating normals" << std::endl;
-   std::vector<Pvl::Vec3f> normals = Pvl::estimateNormals(points);
+   std::vector<Vec3f> normals = estimateNormals(points);
 
    {
        std::ofstream ofs("pc-raw.ply");
-       Pvl::PlyWriter writer(ofs);
+       PlyWriter writer(ofs);
        writer.write(points, normals);
    }
 
    std::cout << "Orienting normals" << std::endl;
-   Pvl::orientNormals(points, normals);
+   orientNormals(points, normals);
    {
        std::ofstream ofs("pc-oriented.ply");
-       Pvl::PlyWriter writer(ofs);
+       PlyWriter writer(ofs);
        writer.write(points, normals);
    }
 */
 /*std::ifstream
 ifs("/home/pavel/projects/random/mesh/data/bunny/reconstruction/bun_zipper.ply");
-Pvl::PlyReader reader(ifs);
+PlyReader reader(ifs);
 auto mesh = reader.readMesh();*/
 
 /* for (int i = 0; i < 200; ++i) {
      std::cout << "iter = " << i << std::endl;
-     Pvl::laplacianSmoothing<Pvl::ParallelTag>(mesh);
+     laplacianSmoothing<ParallelTag>(mesh);
  }
  std::ofstream ofs("bunny.ply");
- Pvl::PlyWriter writer(ofs);
+ PlyWriter writer(ofs);
  writer << mesh;
 
  std::ofstream ofs2("boundary.ply");
- Pvl::PlyWriter boun(ofs2);
+ PlyWriter boun(ofs2);
  for (auto vh : mesh.vertexRange()) {
      if (mesh.boundary(vh)) {
          boun << mesh.points[vh];
      }
  }*/
 
-/*std::vector<Pvl::Vec3f> normals = Pvl::estimateNormals(mesh.points);
+/*std::vector<Vec3f> normals = estimateNormals(mesh.points);
 
 {
     std::ofstream ofs("bunny-raw.ply");
-    Pvl::PlyWriter writer(ofs);
+    PlyWriter writer(ofs);
     writer.write(mesh.points, normals);
 }
 
-Pvl::orientNormals(mesh.points, normals);
+orientNormals(mesh.points, normals);
 {
     std::ofstream ofs("bunny-oriented.ply");
-    Pvl::PlyWriter writer(ofs);
+    PlyWriter writer(ofs);
     writer.write(mesh.points, normals);
 */