changed default index type to uint32_t

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

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) {

Optional.hpp

@@ -1,6 +1,7 @@
 #pragma once
 #include "Assert.hpp"
 #include <type_traits>
+#include <utility>
 
 namespace Pvl {
 

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;
     }

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;
             }

QuadricDecimator.hpp

@@ -28,17 +28,19 @@ 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];
             plane[1] = n[1];
             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;
 

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]);
         });
 }

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;
             }
         }

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

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) {