Add enhanced rsr reconstruction to pygel

src/PyGEL/hmesh_functions.cpp

@@ -10,6 +10,7 @@
 #include <string>
 #include <GEL/HMesh/HMesh.h>
 #include <GEL/HMesh/face_loop.h>
+#include <GEL/HMesh/RSRExperimental.h>
 #include <GEL/Geometry/Graph.h>
 #include <GEL/Geometry/graph_io.h>
 #include <GEL/Geometry/graph_skeletonize.h>
@@ -258,8 +259,67 @@ void rsr_recon(Manifold_ptr m_ptr, double* verts, double* normals, int v_num, in
 
     reconstruct_single(*(reinterpret_cast<Manifold*>(m_ptr)), vertices, norm, 
         isEuclidean, genus, k, r, theta, n);
+}
+
+void rsr_recon_experimental(Manifold_ptr m_ptr, double* verts,
+    double* normals, int v_num, int n_num, bool isEuclidean, int genus,
+    int k, int r, int theta, int n)
+{
+    vector<Vec3d> vertices;
+    vector<Vec3d> norm;
+    vertices.reserve(v_num);
+    norm.reserve(n_num);
+    for (int i = 0; i < v_num; i++) {
+        vertices.emplace_back(verts[i], verts[i + v_num], verts[i + 2 * v_num]);
+    }
+
+    for (int i = 0; i < n_num; i++) {
+        norm.emplace_back(normals[i], normals[i + n_num], normals[i + 2 * n_num]);
+    }
+    RSR::RSROpts opts;
+    opts.dist = (isEuclidean) ? RSR::Distance::Euclidean : RSR::Distance::Tangent;
+    opts.genus = genus;
+    opts.k = k;
+    opts.r = r;
+    opts.theta = theta;
+    opts.n = n;
+
+    Manifold result = RSR::point_cloud_to_mesh(vertices, norm, opts);
+    *reinterpret_cast<Manifold*>(m_ptr) = std::move(result);
+}
+
+void hrsr_recon_experimental(Manifold_ptr m_ptr, double* verts, double* normals, size_t v_num, size_t n_num,
+    int collapse_iters, bool is_euclidean,
+    int genus, int k, int r, int theta, int n,
+    bool skip_reexpansion
+    )
+{
+    vector<Vec3d> vertices;
+    vector<Vec3d> norm;
+    vertices.reserve(v_num);
+    norm.reserve(n_num);
+    for (int i = 0; i < v_num; i++) {
+        vertices.emplace_back(verts[i], verts[i + v_num], verts[i + 2 * v_num]);
+    }
+
+    for (int i = 0; i < n_num; i++) {
+        norm.emplace_back(normals[i], normals[i + n_num], normals[i + 2 * n_num]);
+    }
 
-    return; 
+    RSR::CollapseOpts collapse_opts;
+    collapse_opts.max_iterations = collapse_iters;
+    collapse_opts.distance = (is_euclidean) ? RSR::Distance::Euclidean : RSR::Distance::Tangent;
+    RSR::RSROpts rsr_opts;
+    rsr_opts.genus = genus;
+    rsr_opts.k = k;
+    rsr_opts.r = r;
+    rsr_opts.theta = theta;
+    rsr_opts.n = n;
+    RSR::ReexpandOpts reexpand_opts;
+    reexpand_opts.enabled = !skip_reexpansion;
+
+    Manifold result = RSR::point_cloud_collapse_reexpand(vertices, norm, collapse_opts, rsr_opts, reexpand_opts);
+    *reinterpret_cast<Manifold*>(m_ptr) = std::move(result);
 }
 
 using IntVector = vector<size_t>;

src/PyGEL/hmesh_functions.h

@@ -118,6 +118,16 @@ extern "C" {
         double* normals, int v_num, int n_num, bool isEuclidean = false, int genus = 0,
         int k = 70, int r = 20, int theta = 60, int n = 50);
 
+    DLLEXPORT void rsr_recon_experimental(Manifold_ptr m_ptr, double* verts,
+        double* normals, int v_num, int n_num, bool isEuclidean = false, int genus = 0,
+        int k = 70, int r = 20, int theta = 60, int n = 50);
+
+    DLLEXPORT void hrsr_recon_experimental(Manifold_ptr m_ptr, double* verts, double* normals, size_t v_num, size_t n_num,
+    int collapse_iters = 4, bool is_euclidean = false,
+    int genus = -1, int k = 30, int r = 20, int theta = 60, int n = 50,
+    bool skip_reexpansion = false
+    );
+
     DLLEXPORT void extrude_faces(Manifold_ptr _m_ptr, int* faces, int no_faces, IntVector_ptr _fidx_ptr);
 
     DLLEXPORT void kill_face_loop(Manifold_ptr _m_ptr);

src/PyGEL/pygel3d/__init__.py

@@ -30,7 +30,7 @@
 PyGEL is based on the C++ GEL library and provides a Python interface for most but not
 all of the functionality of GEL. 
 """
-__all__ = ["hmesh", "graph", "gl_display", "jupyter_display", "spatial"]
+__all__ = ["hmesh", "graph", "gl_display", "jupyter_display", "spatial", "experimental"]
 
 import os
 from sys import platform, prefix
@@ -205,6 +205,8 @@ def _get_lib_name():
 lib_py_gel.ply_load.argtypes = (ct.c_char_p, ct.c_void_p)
 lib_py_gel.x3d_load.argtypes = (ct.c_char_p, ct.c_void_p)
 lib_py_gel.rsr_recon.argtypes = (ct.c_void_p, ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ct.c_int, ct.c_int, ct.c_bool, ct.c_int, ct.c_int, ct.c_int, ct.c_int, ct.c_int)
+lib_py_gel.rsr_recon_experimental.argtypes = (ct.c_void_p, ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ct.c_int, ct.c_int, ct.c_bool, ct.c_int, ct.c_int, ct.c_int, ct.c_int, ct.c_int)
+lib_py_gel.hrsr_recon_experimental.argtypes = (ct.c_void_p, ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ndpointer(ndim=2, dtype=ct.c_double,flags='F'), ct.c_size_t, ct.c_size_t, ct.c_int, ct.c_bool, ct.c_int, ct.c_int, ct.c_int, ct.c_int, ct.c_int, ct.c_bool)
 lib_py_gel.remove_caps.argtypes = (ct.c_void_p, ct.c_float)
 lib_py_gel.remove_needles.argtypes = (ct.c_void_p, ct.c_float, ct.c_bool)
 lib_py_gel.close_holes.argtypes = (ct.c_void_p,ct.c_int)

src/PyGEL/pygel3d/experimental/__init__.py

@@ -0,0 +1 @@
+__all__ = ["hmesh"]

src/PyGEL/pygel3d/experimental/hmesh.py

@@ -0,0 +1,58 @@
+import numpy as np
+import ctypes as ct
+from numpy.typing import ArrayLike
+from pygel3d.hmesh import Manifold
+from pygel3d import lib_py_gel
+
+def rsr_recon(verts: ArrayLike,
+              normals: ArrayLike=None,
+              use_Euclidean_distance: bool=False,
+              genus: int=-1,
+              k: int=70,
+              r: float=20,
+              theta: float=60,
+              n: int=50) -> Manifold:
+    """ RsR Reconstruction. The first argument, verts, is the point cloud. The next argument,
+        normals, are the normals associated with the vertices or empty list (default) if normals
+        need to be estimated during reconstruction. use_Euclidean_distance should be true if we
+        can use the Euclidean rather than projected distance. Set to true only for noise free
+        point clouds. genus is used to constrain the genus of the reconstructed object. genus
+        defaults to -1, meaning unknown genus. k is the number of nearest neighbors for each point,
+        r is the maximum distance to farthest neighbor measured in multiples of average distance,
+        theta is the threshold on angles between normals: two points are only connected if the angle
+        between their normals is less than theta. Finally, n is the threshold on the distance between
+        vertices that are connected by handle edges (check paper). For large n, it is harder for
+        the algorithm to add handles. """
+    m = Manifold()
+    verts_data = np.asarray(verts, dtype=ct.c_double, order='F')
+    n_verts = len(verts)
+    n_normal = 0 if normals is None else len(normals)
+    if(n_normal==0):
+        normals = [[]]
+    normal_data = np.asarray(normals, dtype=ct.c_double, order='F')
+
+    lib_py_gel.rsr_recon_experimental(m.obj, verts_data, normal_data, n_verts, n_normal,
+                         use_Euclidean_distance, genus, k, r, theta, n)
+    return m
+
+def hrsr_recon(verts: ArrayLike,
+              normals: ArrayLike=None,
+              collapse_iters = 4,
+              use_Euclidean_distance: bool=False,
+              genus: int=-1,
+              k: int=70,
+              r: float=20,
+              theta: float=60,
+              n: int=50,
+              skip_reexpansion = False) -> Manifold:
+    m = Manifold()
+    verts_data = np.asarray(verts, dtype=ct.c_double, order='F')
+    n_verts = len(verts)
+    n_normal = 0 if normals is None else len(normals)
+    if(n_normal==0):
+        normals = [[]]
+    normal_data = np.asarray(normals, dtype=ct.c_double, order='F')
+
+    lib_py_gel.hrsr_recon_experimental(m.obj, verts_data, normal_data, n_verts, n_normal,
+                                      collapse_iters, use_Euclidean_distance, genus, k, r, theta, n, skip_reexpansion)
+    return m