diff --git a/src/gpytoolbox/__init__.py b/src/gpytoolbox/__init__.py index 082e1d79..1c1c2528 100644 --- a/src/gpytoolbox/__init__.py +++ b/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 diff --git a/src/gpytoolbox/poisson_surface_reconstruction.py b/src/gpytoolbox/stochastic_poisson_surface_reconstruction.py similarity index 94% rename from src/gpytoolbox/poisson_surface_reconstruction.py rename to src/gpytoolbox/stochastic_poisson_surface_reconstruction.py index e60d8a73..4ea20d71 100644 --- a/src/gpytoolbox/poisson_surface_reconstruction.py +++ b/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 diff --git a/test/test_point_cloud_to_mesh.py b/test/test_point_cloud_to_mesh.py index 368b8211..0f33b6b8 100644 --- a/test/test_point_cloud_to_mesh.py +++ b/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__': diff --git a/test/test_poisson_surface_reconstruction.py b/test/test_stochastic_poisson_surface_reconstruction.py similarity index 93% rename from test/test_poisson_surface_reconstruction.py rename to test/test_stochastic_poisson_surface_reconstruction.py index 8da38067..82034b7b 100644 --- a/test/test_poisson_surface_reconstruction.py +++ b/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)