Poisson surface reconstruction volume

I have surfaces that are represented by noisy point clouds that I am trying to extract an accurate volume from.

The actual surface can be assumed to be somewhere within the noisy point cloud surface (cannot share the data but can create an example if needed). A Poisson surface reconstruction is a decent way to get a realistic surface from these point clouds.

I am able to do a Poisson surface reconstruction and get the volume contained (if mesh is watertight) in an external program like MeshLab. This is the volume I am hoping to achieve, without having to leave MATLAB since I need to calculate many volumes in my workflow.

Things I've tried:

• Convex hull - easy and fast to calculate but leads to an overestimate since it is sensitive to outlier points.
• Alpha shape - essentially gets inner volume. Works ok unless there are holes in the point cloud (which is often), in which case it is wildly inaccurate.
• Converting a point cloud to a mesh (using pc2surfacemesh) then back to a point cloud (using mesh2pc) and measuring volume via convex hull - this is the most accurate approach so far, but still leads to volume overestimation compared to using MeshLab.

In summary, is there any way to get a volume directly from a Poisson surface reconstruction in MATLAB?

Answer 1

[April A]

My solution is going to be using PyMeshLab to do the Poisson surface reconstruction and get the volume in a Python script, then calling that script from MATLAB.

Answer 2

[Prathamesh]

Hi April A,

I understand that you want to get a volume directly from a Poisson surface reconstruction in MATLAB. You have surfaces that are represented by noisy point clouds that you are trying to extract an accurate volume from.

First of all, directly implementing Poisson surface reconstruction in MATLAB is challenging because MATLAB does not have built-in support for this specific algorithm. However, you can achieve a similar result through a combination of MATLAB functions and by approximating the process.

So, we can do this by creating a voxel grid and then using “ isosurface ” function of MATLAB. You can create a voxel grid and use “isosurface” to generate a mesh. This method approximates the surface but may not be as precise as Poisson reconstruction.

Below is an example of code use to demonstrate above idea,

% Generating sample point cloud data resembling a noisy sphere
numPoints = 1000;
theta = 2 * pi * rand(numPoints, 1);
phi = acos(2 * rand(numPoints, 1) - 1);
r = 1 + 0.05 * randn(numPoints, 1); % Add some noise
x = r .* sin(phi) .* cos(theta);
y = r .* sin(phi) .* sin(theta);
z = r .* cos(phi);
pointCloudData = [x, y, z];

% Defining the voxel grid resolution
voxelSize = 0.05;
minBound = min(pointCloudData) - voxelSize;
maxBound = max(pointCloudData) + voxelSize;
gridSize = ceil((maxBound - minBound) / voxelSize);

% Calculate voxel indices
indices = floor((pointCloudData - minBound) / voxelSize) + 1;

% Create a 3D histogram (voxel grid) using accumarray
voxelGrid = accumarray(indices, 1, gridSize);

% Generate an isosurface from the voxel grid
[X, Y, Z] = ndgrid(...
    minBound(1):voxelSize:maxBound(1), ...
    minBound(2):voxelSize:maxBound(2), ...
    minBound(3):voxelSize:maxBound(3));
isoValue = 0.5; % Adjust this value based on your data
surfaceMesh = isosurface(X, Y, Z, voxelGrid, isoValue);

% Plot the resulting mesh
figure;
patch(surfaceMesh, 'FaceColor', [0.8, 0.8, 1.0], 'EdgeColor', 'none');
camlight; lighting phong;
xlabel('X'); ylabel('Y'); zlabel('Z');
title('Reconstructed Surface Mesh');
axis equal;

% Calculate the volume of the mesh
volume = meshVolume(surfaceMesh.vertices, surfaceMesh.faces);
fprintf('Estimated Volume: %.3f\n', volume);

% Helper function to calculate volume of a mesh
function vol = meshVolume(vertices, faces)
vol = 0;
for i = 1:size(faces, 1)
    v0 = vertices(faces(i, 1), :);
    v1 = vertices(faces(i, 2), :);
    v2 = vertices(faces(i, 3), :);
    tetraVolume = dot(v0, cross(v1, v2)) / 6;
    vol = vol + tetraVolume;
end
vol = abs(vol);
end

I have attached a screenshot of the output.

[April A]

Hi Prathamesh,

Thanks for your answer and helpful example.

Can you verify that your code is giving the correct answer? In your example, you generated a sphere with radius = 1, which should have a volume (4/3 * pi * r^3) of about 4.19.

If your screenshot and when I run the code myself, the volume is estimated to be 0.028.

Additionally, the surface looks quite patchy leading me to think the volume calculated is not the volume interior to the sphere.

[Prathamesh]

Hi @April A can you please let me know what answer are you expecting after running the above code I provided? Because my code is calulation the volume of the patchy cloudy surface.

[April A]

@Prathamesh I would expect the volume of a sphere with radius 1 to be ~4.19.

I am interested in calculating the volume enclosed by the surface, not the volume of the surface itself.

Answer 3

[Prathamesh]

@April A Please let me know if this works

import pymeshlab
import numpy as np

# Create a mesh from the noisy sphere data
num_points = 1000
theta = 2 * np.pi * np.random.rand(num_points)
phi = np.arccos(2 * np.random.rand(num_points) - 1)
r = 1 + 0.05 * np.random.randn(num_points)
x = r * np.sin(phi) * np.cos(theta)
y = r * np.sin(phi) * np.sin(theta)
z = r * np.cos(phi)
point_cloud_data = np.column_stack((x, y, z))

# Initialize a new mesh from the point cloud
ms = pymeshlab.MeshSet()
ms.add_mesh(pymeshlab.Mesh(vertex_matrix=point_cloud_data))

# Apply Poisson surface reconstruction
ms.apply_filter('compute_normals_for_point_sets')
ms.apply_filter('surface_reconstruction_screened_poisson', preclean=True)

# Calculate the volume of the reconstructed mesh
volume = ms.compute_geometric_measures().mesh_volume
print(f'Estimated Volume: {volume:.3f}')

To call this Python script from MATLAB, you can use the system command or the py module if you have Python integrated into MATLAB:

system('python poisson_reconstruction.py');

Or using the py module:

py.poisson_reconstruction.main();

[April A]

Correct, this is what I'm doing, more or less.

pyrunfile works well in MATLAB.