Update README so examples don't have any errors. Add a few tests. Fix…

… some small bugs.

CMakeLists.txt

@@ -145,7 +145,7 @@ npe_add_module(_pcu_internal
   ${CMAKE_CURRENT_SOURCE_DIR}/src/curvature.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/sparse_voxel_grid.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/marching_cubes.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/src/decimate_mesh.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/mesh_decimate.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/remove_unreferenced_mesh_vertices.cpp
   EXTRA_MODULE_FUNCTIONS
   hack_extra_bindings

README.md

@@ -92,6 +92,7 @@ The following dependencies are required to install with `pip`:
 - [Making a mesh watertight](#making-a-mesh-watertight)
 - [Ray/Mesh intersection](#ray-mesh-intersection)
 - [Ray/Surfel intersection](#ray-surfel-intersection)
+- [Computing curvature on a mesh](#computing-curvature-on-a-mesh)
 
 
 ### Loading meshes and point clouds
@@ -162,6 +163,8 @@ mesh = pcu.TriangleMesh("path/to/mesh")
 For meshes and point clouds with more complex attributes, use `save_triangle_mesh` which accepts a whole host of named
 arguments which control the attributes to save.
 ```python
+import point_cloud_utils as pcu
+
 # save_triangle_mesh accepts a path to save to (The type of mesh  saved is determined by the file extesion),
 # an array of mesh vertices of shape [V, 3], and optional arguments specifying faces, per-mesh attributes,
 # per-face attributes and per-wedge attributes:
@@ -231,7 +234,6 @@ v, f, n, c = pcu.save_mesh_vfnc("path/to/mesh", v, f, n, c)
 Generate 10000 samples on a mesh with poisson disk samples
 ```python
 import point_cloud_utils as pcu
-import numpy as np
 
 # v is a nv by 3 NumPy array of vertices
 # f is an nf by 3 NumPy array of face indexes into v
@@ -243,8 +245,8 @@ v, f, n = pcu.load_mesh_vfn("my_model.ply")
 f_i, bc = pcu.sample_mesh_poisson_disk(v, f, n, 10000)
 
 # Use the face indices and barycentric coordinate to compute sample positions and normals
-v_poisson = pcu.interpolate_barycentric_coords(f, fi, bc, v)
-n_poisson = pcu.interpolate_barycentric_coords(f, fi, bc, n)
+v_poisson = pcu.interpolate_barycentric_coords(f, f_i, bc, v)
+n_poisson = pcu.interpolate_barycentric_coords(f, f_i, bc, n)
 ```
 
 Generate blue noise samples on a mesh separated by approximately 0.01 times the bounding box diagonal
@@ -266,8 +268,8 @@ bbox_diag = np.linalg.norm(bbox)
 f_i, bc = pcu.sample_mesh_poisson_disk(v, f, n, 10000)
 
 # Use the face indices and barycentric coordinate to compute sample positions and normals
-v_sampled = pcu.interpolate_barycentric_coords(f, fi, bc, v)
-n_sampled = pcu.interpolate_barycentric_coords(f, fi, bc, n)
+v_sampled = pcu.interpolate_barycentric_coords(f, f_i, bc, v)
+n_sampled = pcu.interpolate_barycentric_coords(f, f_i, bc, n)
 ```
 
 ### Generate random samples on a mesh
@@ -281,12 +283,12 @@ import numpy as np
 v, f, n = pcu.load_mesh_vfn("my_model.ply")
 
 # Generate random samples on the mesh (v, f, n)
-# f_idx are the face indices of each sample and bc are barycentric coordinates of the sample within a face
-f_idx, bc = pcu.sample_mesh_random(v, f, num_samples=v.shape[0] * 40)
+# f_i are the face indices of each sample and bc are barycentric coordinates of the sample within a face
+f_i, bc = pcu.sample_mesh_random(v, f, num_samples=v.shape[0] * 40)
 
 # Use the face indices and barycentric coordinate to compute sample positions and normals
-v_sampled = pcu.interpolate_barycentric_coords(f, fi, bc, v)
-n_sampled = pcu.interpolate_barycentric_coords(f, fi, bc, n)
+v_sampled = pcu.interpolate_barycentric_coords(f, f_i, bc, v)
+n_sampled = pcu.interpolate_barycentric_coords(f, f_i, bc, n)
 ```
 
 ### Downsample a point cloud to have a blue noise distribution
@@ -395,6 +397,7 @@ v_sampled, n_sampled, c_sampled = pcu.downsample_point_cloud_voxel_grid(sizeof_v
 ### Compute closest points on a mesh
 ```python
 import point_cloud_utils as pcu
+import numpy as np
 
 # v is a nv by 3 NumPy array of vertices
 v, f = pcu.load_mesh_vf("my_model.ply")
@@ -499,14 +502,16 @@ import point_cloud_utils as pcu
 import numpy as np
 
 # Generate two random point sets
-a = np.random.rand(1000, 3)
-b = np.random.rand(500, 3)
+pts_a = np.random.rand(1000, 3)
+pts_b = np.random.rand(500, 3)
+
+k = 10
 
-# dists_a_to_b is of shape (a.shape[0],) and contains the shortest squared distance
-# between each point in a and the points in b
-# corrs_a_to_b is of shape (a.shape[0],) and contains the index into b of the
-# closest point for each point in a
-dists_a_to_b, corrs_a_to_b = pcu.shortest_distance_pairs(a, b)
+# dists_a_to_b is of shape (pts_a.shape[0], k) and contains the (sorted) distances 
+# to the k nearest points in pts_b
+# corrs_a_to_b is of shape (a.shape[0], k) and contains the index into pts_b of the 
+# k closest points for each point in pts_a 
+dists_a_to_b, corrs_a_to_b = pcu.k_nearest_neighbors(pts_a, pts_b, k)
 ```
 
 ### Generating point samples in the square and cube with Lloyd relaxation
@@ -530,6 +535,7 @@ samples_3d = pcu.lloyd_3d(100)
 ### Compute shortest signed distances to a triangle mesh with [fast winding numbers](https://www.dgp.toronto.edu/projects/fast-winding-numbers/)
 ```python
 import point_cloud_utils as pcu
+import numpy as np
 
 # v is a nv by 3 NumPy array of vertices
 # f is an nf by 3 NumPy array of face indexes into v
@@ -557,7 +563,7 @@ p, n = pcu.load_mesh_vn("my_pcloud.ply")
 # Treat any points closer than 1e-7 apart as the same point
 # idx_i is an array of indices such that p_dedup = p[idx_i]
 # idx_j is an array of indices such that p = p_dedup[idx_j]
-p_dedup, idx_i, idx_j  = deduplicate_point_cloud(p, 1e-7)
+p_dedup, idx_i, idx_j  = pcu.deduplicate_point_cloud(p, 1e-7)
 
 # Use idx_i to deduplicate the normals
 n_dedup = n[idx_i]
@@ -661,7 +667,7 @@ v_watertight, f_watertight = pcu.make_mesh_watertight(v, f, resolution=resolutio
 ```
 
 
-# Ray/Mesh Intersection
+### Ray/Mesh Intersection
 ```python
 import point_cloud_utils as pcu
 import numpy as np
@@ -675,7 +681,7 @@ v, f, c = pcu.load_mesh_vfc("my_model.ply")
 uv = np.stack([a.ravel() for a in np.mgrid[-1:1:128j, -1.:1.:128j]], axis=-1)
 ray_d = np.concatenate([uv, np.ones([uv.shape[0], 1])], axis=-1)
 ray_d = ray_d / np.linalg.norm(ray_d, axis=-1, keepdims=True)
-ray_o = np.array([[2.5, 0, -55.0] for _ in range(d.shape[0])])
+ray_o = np.array([[2.5, 0, -55.0] for _ in range(ray_d.shape[0])])
 
 # Intersect rays with geometry
 intersector = pcu.RayMeshIntersector(v, f)
@@ -691,7 +697,7 @@ hit_pos = pcu.interpolate_barycentric_coords(f, fid[hit_mask], bc[hit_mask], v)
 hit_clr = pcu.interpolate_barycentric_coords(f, fid[hit_mask], bc[hit_mask], c)
 ```
 
-# Ray/Surfel Intersection
+### Ray/Surfel Intersection
 ```python
 import point_cloud_utils as pcu
 import numpy as np
@@ -704,7 +710,7 @@ v, n = pcu.load_mesh_vn("my_model.ply")
 uv = np.stack([a.ravel() for a in np.mgrid[-1:1:128j, -1.:1.:128j]], axis=-1)
 ray_d = np.concatenate([uv, np.ones([uv.shape[0], 1])], axis=-1)
 ray_d = ray_d / np.linalg.norm(ray_d, axis=-1, keepdims=True)
-ray_o = np.array([[2.5, 0, -55.0] for _ in range(d.shape[0])])
+ray_o = np.array([[2.5, 0, -55.0] for _ in range(ray_d.shape[0])])
 
 # Intersect rays with surfels with fixed radius 0.55
 intersector = pcu.RaySurfelIntersector(v, n, r=0.55)
@@ -717,3 +723,27 @@ pid, t = intersector.intersect_rays(ray_o, ray_d)
 hit_mask = pid >= 0
 intersected_points = v[pid[hit_mask]]
 ```
+
+### Computing curvature on a mesh
+```python
+import point_cloud_utils as pcu
+
+# v is a #v by 3 NumPy array of vertices
+# f is an #f by 3 NumPy array of face indexes into v
+v, f = pcu.load_mesh_vfc("my_model.ply")
+
+# Compute principal min/max curvature magnitudes (k1, k2) and directions (d1, d2) 
+# using the one ring of each vertex
+k1, k2, d1, d2 = pcu.mesh_principal_curvatures(v, f)
+
+# Compute principal min/max curvature magnitudes (k1, k2) and directions (d1, d2) 
+# using a radius. This method is much more robust but requires tuning the radius
+k1, k2, d1, d2 = pcu.mesh_principal_curvatures(v, f, r=0.1)
+
+# Compute Mean (kh) and Gaussian (kg) curvatures using the one ring of each vertex
+kh, kg = pcu.mesh_mean_and_gaussian_curvatures(v, f)
+
+# Compute Mean (kh) and Gaussian (kg) curvatures using using a radius.
+# This method is much more robust but requires tuning the radius
+kh, kg = pcu.mesh_mean_and_gaussian_curvatures(v, f, r=0.1)
+```

point_cloud_utils/__init__.py

@@ -5,7 +5,7 @@
     lloyd_2d, lloyd_3d, voronoi_centroids_unit_cube, sample_mesh_lloyd, \
     deduplicate_point_cloud, deduplicate_mesh_vertices, signed_distance_to_mesh, \
     closest_points_on_mesh, connected_components, ray_mesh_intersection, laplacian_smooth_mesh, \
-    make_mesh_watertight, mesh_principal_curvatures, mesh_mean_and_gaussian_curvatures, \
+    make_mesh_watertight, mesh_principal_curvatures, \
     morton_add, morton_subtract, \
     sparse_voxel_grid_boundary, marching_cubes_sparse_voxel_grid, decimate_triangle_mesh, \
     remove_unreferenced_mesh_vertices
@@ -18,6 +18,29 @@
 from ._ray_point_cloud_intersector import ray_surfel_intersection, surfel_geometry, RaySurfelIntersector
 
 
+def mesh_mean_and_gaussian_curvatures(v, f, r=-1.0):
+    """
+    Estimate mean and Gaussian curvatures for a mesh
+
+    Parameters
+    ----------
+    v : #v by 3 Matrix of mesh vertex 3D positions
+    f : #f by 3 Matrix of face (triangle) indices
+    r : optional floating point radius of neighborhood to consider when estimating curvature
+        If set to a positive value, will use a more robust curvature estimation method (but may require some tuning)
+
+    Returns
+    -------
+    A tuple (kh, kg) where:
+    kh is an array of shape (#v,) of per-vertex mean curvatures
+    kg is an array of shape (#v,) of per-vertex Gaussian curvatures
+    :return:
+    """
+    k1, k2, _, _ = mesh_principal_curvatures(v, f, r)
+
+    return 0.5 * (k1 + k2), (k1 * k2)
+
+
 def hausdorff_distance(x, y, return_index=False, squared_distances=False, max_points_per_leaf=10):
     """
     Compute the Hausdorff distance between x and y

src/curvature.cpp

@@ -81,44 +81,44 @@ npe_begin_code()
 npe_end_code()
 
 
-const char* mesh_mean_and_gaussian_curvatures_doc = R"Qu8mg5v7(
-Estimate mean and Gaussian curvatures for a mesh
-
-Parameters
-----------
-v : #v by 3 Matrix of mesh vertex 3D positions
-f : #f by 3 Matrix of face (triangle) indices
-
-Returns
--------
-A tuple (kg, kh, d1, d2) where:
-  kg is an array of shape (#v,) of per-vertex Gaussian curvatures
-  kh is an array of shape (#v,) of per-vertex mean curvatures
-)Qu8mg5v7";
-npe_function(mesh_mean_and_gaussian_curvatures)
-npe_arg(v, dense_float, dense_double)
-npe_arg(f, dense_int, dense_long, dense_longlong)
-npe_doc(mesh_mean_and_gaussian_curvatures_doc)
-npe_begin_code()
-{
-    VCGMesh m;
-    vcg_mesh_from_vf(v, f, m);
-	tri::UpdateTopology<VCGMesh>::FaceFace(m);
-	tri::UpdateTopology<VCGMesh>::VertexFace(m);
-    tri::UpdateCurvature<VCGMesh>::MeanAndGaussian(m);
-
-    npe_Matrix_v kg(m.vn, 1), kh(m.vn, 1);
-
-    int vcount = 0;
-    for (VCGMesh::VertexIterator vit = m.vert.begin(); vit != m.vert.end(); vit++) {
-        kg(vcount, 0) = vit->cKg();
-        kh(vcount, 0) = vit->cKh();
-        vcount += 1;
-    }
-
-    return std::make_tuple(npe::move(kg), npe::move(kh));
-}
-npe_end_code()
+//const char* mesh_mean_and_gaussian_curvatures_doc = R"Qu8mg5v7(
+//Estimate mean and Gaussian curvatures for a mesh
+//
+//Parameters
+//----------
+//v : #v by 3 Matrix of mesh vertex 3D positions
+//f : #f by 3 Matrix of face (triangle) indices
+//
+//Returns
+//-------
+//A tuple (kg, kh, d1, d2) where:
+//  kg is an array of shape (#v,) of per-vertex Gaussian curvatures
+//  kh is an array of shape (#v,) of per-vertex mean curvatures
+//)Qu8mg5v7";
+//npe_function(mesh_mean_and_gaussian_curvatures)
+//npe_arg(v, dense_float, dense_double)
+//npe_arg(f, dense_int, dense_long, dense_longlong)
+//npe_doc(mesh_mean_and_gaussian_curvatures_doc)
+//npe_begin_code()
+//{
+//    VCGMesh m;
+//    vcg_mesh_from_vf(v, f, m);
+//	tri::UpdateTopology<VCGMesh>::FaceFace(m);
+//	tri::UpdateTopology<VCGMesh>::VertexFace(m);
+//    tri::UpdateCurvature<VCGMesh>::MeanAndGaussian(m);
+//
+//    npe_Matrix_v kg(m.vn, 1), kh(m.vn, 1);
+//
+//    int vcount = 0;
+//    for (VCGMesh::VertexIterator vit = m.vert.begin(); vit != m.vert.end(); vit++) {
+//        kg(vcount, 0) = vit->cKg();
+//        kh(vcount, 0) = vit->cKh();
+//        vcount += 1;
+//    }
+//
+//    return std::make_tuple(npe::move(kg), npe::move(kh));
+//}
+//npe_end_code()
 
 
 //const char* pointcloud_apss_curvature_doc = R"Qu8mg5v7(

src/decimate_mesh.cpp → src/mesh_decimate.cpp

@@ -1,6 +1,5 @@
 #include <npe.h>
 #include <igl/decimate.h>
-#include <igl/remove_unreferenced.h>
 
 
 #include "common/common.h"
@@ -63,7 +62,7 @@ npe_begin_code()
         out_face_correspondences = out_f_corr_copy.template cast<npe_Scalar_f>();
         out_vertex_correspondences = out_v_corr_copy.template cast<npe_Scalar_f>();
 
-        return std::make_tuple(npe::move(out_v), npe::move(out_f), npe::move(out_face_correspondences), npe::move(out_vertex_correspondences));
+        return std::make_tuple(npe::move(out_v), npe::move(out_f), npe::move(out_vertex_correspondences), npe::move(out_face_correspondences));
     } else {
         throw pybind11::value_error("Invalid decimation heuristic, must be 'shortest_edge'. Others will be implemented"
                                     " in future versions of Point Cloud Utils.");

tests/test_examples.py

@@ -590,6 +590,28 @@ def test_ply_load_save(self):
         v, f = pcu.load_mesh_vf(os.path.join(self.test_path, "duplicated_pcloud.ply"))
         pcu.save_mesh_vf("test.ply", v, f)
 
+    def test_mesh_curvature(self):
+        import point_cloud_utils as pcu
+        import numpy as np
+        v, f = pcu.load_mesh_vf(os.path.join(self.test_path, "bunny.ply"))
+        k1, k2, d1, d2 = pcu.mesh_principal_curvatures(v, f)
+        kh, kg = pcu.mesh_mean_and_gaussian_curvatures(v, f)
+
+        self.assertTrue(np.allclose(0.5 * (k1 + k2), kh))
+        self.assertTrue(np.allclose((k1 * k2), kg))
+
+    def test_mesh_decimation(self):
+        import point_cloud_utils as pcu
+        v, f = pcu.load_mesh_vf(os.path.join(self.test_path, "bunny.ply"))
+        target_num_faces = f.shape[0] // 4
+
+        v_decimate, f_decimate, v_correspondence, f_correspondence = pcu.decimate_triangle_mesh(v, f, target_num_faces)
+        birth_v, birth_f = v[v_correspondence], f[f_correspondence]
+
+        self.assertLess(v_correspondence.max(), v.shape[0])
+        self.assertLess(f_correspondence.max(), f.shape[0])
+        self.assertLessEqual(f_decimate.shape[0], target_num_faces)
+
 
 if __name__ == '__main__':
     unittest.main()