stochastic psr rework, fixed test

src/gpytoolbox/__init__.py

@@ -13,7 +13,7 @@
 from .signed_distance_polygon import signed_distance_polygon
 from .metropolis_hastings import metropolis_hastings
 from .ray_polyline_intersect import ray_polyline_intersect
-from .poisson_surface_reconstruction import poisson_surface_reconstruction
+from .stochastic_poisson_surface_reconstruction import stochastic_poisson_surface_reconstruction
 from .fd_interpolate import fd_interpolate
 from .fd_grad import fd_grad
 from .fd_partial_derivative import fd_partial_derivative

src/gpytoolbox/poisson_surface_reconstruction.py → src/gpytoolbox/stochastic_poisson_surface_reconstruction.py

@@ -9,9 +9,9 @@
 from .grid_neighbors import grid_neighbors
 from .grid_laplacian_eigenfunctions import grid_laplacian_eigenfunctions
 
-def poisson_surface_reconstruction(P,N,gs=None,h=None,corner=None,stochastic=False,sigma_n=0.0,sigma=0.05,solve_subspace_dim=0,verbose=False,prior_fun=None):
+def stochastic_poisson_surface_reconstruction(P,N,gs=None,h=None,corner=None,output_variance=False,sigma_n=0.0,sigma=0.05,solve_subspace_dim=0,verbose=False,prior_fun=None):
     """
-    Runs Poisson Surface Reconstruction from a set of points and normals to output a scalar field that takes negative values inside the surface and positive values outside the surface.
+    Runs Stochastic Poisson Surface Reconstruction from a set of points and normals to output a scalar field that takes negative values inside the surface and positive values outside the surface.
     
     Parameters
     ----------
@@ -25,8 +25,8 @@ def poisson_surface_reconstruction(P,N,gs=None,h=None,corner=None,stochastic=Fal
         Grid spacing in each dimension
     corner : (dim,) numpy array
         Coordinates of the lower left corner of the grid
-    stochastic : bool, optional (default False)
-        Whether to use "Stochastic Poisson Surface Reconstruction" to output a mean and variance scalar field instead of just one scalar field
+    output_variance : bool, optional (default False)
+        Whether to use to output a mean *and* variance scalar field instead of just the mean scalar field
     sigma_n : float, optional (default 0.0)
         Noise level in the normals
     sigma : float, optional (default 0.05)
@@ -43,17 +43,14 @@ def poisson_surface_reconstruction(P,N,gs=None,h=None,corner=None,stochastic=Fal
     scalar_mean : (gs[0],gs[1],...,gs[dim-1]) numpy array
         Mean of the reconstructed scalar field
     scalar_variance : (gs[0],gs[1],...,gs[dim-1]) numpy array
-        Variance of the reconstructed scalar field. This will only be part of the output if stochastic=True.
+        Variance of the reconstructed scalar field. This will only be part of the output if output_variance=True.
     grid_vertices : list of (gs[0],gs[1],...,gs[dim-1],dim) numpy arrays
         Grid vertices (each element in the list is one dimension), as in the output of np.meshgrid
 
     
     Notes
     -----
-    If you are looking for an efficient implementation of screened Poisson
-    surface reconstruction, please use the function `point_cloud_to_mesh` - this
-    function is a highly customizable implementation of PSR that is slower and
-    operates on a uniform grid.
+    This algorithm implements "Stochastic Poisson Surface Reconstruction" as introduced by Sellán and Jacobson in 2022. If you are only looking to reconstruct a mesh from a point cloud, use the traditional "(Screened) Poisson Surface Reconstruction" by Kazhdan et al. implemented in `point_cloud_to_mesh`.
 
     See [this jupyter notebook](https://colab.research.google.com/drive/1DOXbDmqzIygxoQ6LeX0Ewq7HP4201mtr?usp=sharing) for a full tutorial on how to use this function.
 
@@ -272,13 +269,13 @@ def kernel_fun(X,Y):
             # print(prior_fun(staggered_grid_vertices)[:,dd][:,None].shape)
             # print(N[:,dd][:,None].shape)
             means.append(prior_fun(staggered_grid_vertices)[:,dd][:,None] + k2.T @ K3_inv @ (N[:,dd][:,None] - prior_fun(P)[:,dd][:,None]))
-        if stochastic:
+        if output_variance:
             covs.append(k1 - k2.T @ K3_inv @ k2)
 
 
     # Concatenate the mean and covariance matrices
     vector_mean = np.concatenate(means)
-    if stochastic:
+    if output_variance:
         vector_cov = sigma*sigma*block_diag(covs)
 
     if verbose:
@@ -310,7 +307,7 @@ def kernel_fun(X,Y):
         t10 = time.time()
 
 
-    if stochastic:
+    if output_variance:
         # In this case, we want to compute the variances of the scalar field        
         cov_divergence = G.T @ vector_cov @ G
 

test/test_point_cloud_to_mesh.py

@@ -30,7 +30,7 @@ def test_variety_of_meshes(self):
             for mesh in meshes:
                 V0,F0 = gpy.read_mesh("test/unit_tests_data/" + mesh)
                 V0 = gpy.normalize_points(V0)
-                P,I,u = gpy.random_points_on_mesh(V0, F0, 6*V0.shape[0], rng=rng,
+                P,I,u = gpy.random_points_on_mesh(V0, F0, 20*V0.shape[0], rng=rng,
                     return_indices=True)
                 N = gpy.per_face_normals(V0,F0)[I,:]
 
@@ -39,6 +39,7 @@ def test_variety_of_meshes(self):
                     psr_outer_boundary_type=bdry_type)
 
                 h = gpy.approximate_hausdorff_distance(V0,F0,V,F)
+                print(h)
                 self.assertTrue(h<0.01)
 
 if __name__ == '__main__':

test/test_poisson_surface_reconstruction.py → test/test_stochastic_poisson_surface_reconstruction.py

@@ -7,7 +7,7 @@
 # from mpl_toolkits.axes_grid1 import make_axes_locatable
 # import polyscope as ps
 
-class TestPoissonSurfaceReconstruction(unittest.TestCase):
+class TestStochasticPoissonSurfaceReconstruction(unittest.TestCase):
     def test_paper_figure(self):
         poly = gpytoolbox.png2poly("test/unit_tests_data/illustrator.png")[0]
         poly = poly - np.min(poly)
@@ -30,7 +30,7 @@ def test_paper_figure(self):
         gs = np.array([50,50])
         # corner = np.array([-1,-1])
         # h = np.array([0.04,0.04])
-        scalar_mean, scalar_var, grid_vertices = gpytoolbox.poisson_surface_reconstruction(P,N,gs=gs,solve_subspace_dim=0,verbose=False,stochastic=True)
+        scalar_mean, scalar_var, grid_vertices = gpytoolbox.stochastic_poisson_surface_reconstruction(P,N,gs=gs,solve_subspace_dim=0,verbose=False,output_variance=True)
         # corner = P.min(axis=0)
         # h = (P.max(axis=0) - P.min(axis=0))/gs
 
@@ -85,7 +85,7 @@ def test_indicator(self):
         # Normals are the same as positions on a circle
         N = np.concatenate((np.cos(th),np.sin(th)),axis=1)
         gs = np.array([100,100])
-        scalar_mean, scalar_var, grid_vertices = gpytoolbox.poisson_surface_reconstruction(P,N,gs=gs,solve_subspace_dim=1000,verbose=False,stochastic=True)
+        scalar_mean, scalar_var, grid_vertices = gpytoolbox.stochastic_poisson_surface_reconstruction(P,N,gs=gs,solve_subspace_dim=1000,verbose=False,output_variance=True)
         prob_out = 1 - norm.cdf(scalar_mean,0,np.sqrt(scalar_var))
 
         # The first test we can run is generate many points inside the shape
@@ -142,7 +142,7 @@ def test_3d(self):
             # Windows machine in github action can't handle this test
             pass
         else:
-            scalar_mean, scalar_var, grid_vertices = gpytoolbox.poisson_surface_reconstruction(P,N,corner=np.array([-1.1,-1.1,-1.1]),h=np.array([0.05,0.05,0.05]),gs=gs,solve_subspace_dim=3000,stochastic=True,verbose=False)
+            scalar_mean, scalar_var, grid_vertices = gpytoolbox.stochastic_poisson_surface_reconstruction(P,N,corner=np.array([-1.1,-1.1,-1.1]),h=np.array([0.05,0.05,0.05]),gs=gs,solve_subspace_dim=3000,output_variance=True,verbose=False)
             grid_vertices = np.array(grid_vertices).reshape(3, -1,order='F').T
             # Where is the highest variance?
             variance_argsort = np.argsort(scalar_var)