As of now there is no inbuilt MATLAB function for the computation and plotting of Gaussian curvature. However, here is a sample code for computing the Gaussian curvature for each point in the point cloud.
% Generate synthetic point cloud data
num_points = 1000;
theta = 2 * pi * rand(num_points, 1);
phi = pi * rand(num_points, 1);
r = 1; % Radius of the sphere
% Generate points on a sphere
X = r * sin(phi) .* cos(theta);
Y = r * sin(phi) .* sin(theta);
Z = r * cos(phi);
pointCloudData = [X, Y, Z];
% Create a point cloud object
ptCloud = pointCloud(pointCloudData);
% Use alphaShape to create a mesh from the point cloud
alpha = 1.25 * r; % Adjust alpha as needed
shp = alphaShape(pointCloudData, alpha);
% Extract the surface mesh
[tri, pts] = boundaryFacets(shp);
% Calculate Gaussian curvature using the mesh
% K stores the Gaussian curvature computation for each point
K = computeGaussianCurvature(tri, pts);
% Helper function to compute Gaussian curvature
function K = computeGaussianCurvature(tri, pts)
% Initialize Gaussian curvature array
K = zeros(size(pts, 1), 1);
% Calculate angle deficit and area for each vertex
for i = 1:size(pts, 1)
% Find triangles attached to vertex i
attachedTriangles = find(any(tri == i, 2));
% Initialize angle sum and area sum
angleSum = 0;
areaSum = 0;
% Process each attached triangle
for j = 1:length(attachedTriangles)
% Indices of the vertices of the current triangle
triIndices = tri(attachedTriangles(j), :);
% Calculate angles and area of the current triangle
angles = computeAngles(pts(triIndices, :));
area = triangleArea(pts(triIndices, :));
% Sum angles and areas
angleSum = angleSum + angles(triIndices == i);
areaSum = areaSum + area;
end
% Calculate Gaussian curvature
K(i) = (2 * pi - angleSum) / areaSum;
end
end
% Compute angles of a triangle given its vertices
function angles = computeAngles(vertices)
% Calculate vectors for the edges
v1 = vertices(2, :) - vertices(1, :);
v2 = vertices(3, :) - vertices(1, :);
v3 = vertices(3, :) - vertices(2, :);
% Calculate angles using dot product
angles(1) = acos(dot(v1, v2) / (norm(v1) * norm(v2)));
angles(2) = acos(dot(-v1, v3) / (norm(v1) * norm(v3)));
angles(3) = acos(dot(-v2, -v3) / (norm(v2) * norm(v3)));
end
% Calculate the area of a triangle given its vertices
function area = triangleArea(vertices)
% Calculate vectors for the edges
v1 = vertices(2, :) - vertices(1, :);
v2 = vertices(3, :) - vertices(1, :);
% Calculate area using cross product
area = 0.5 * norm(cross(v1, v2));
end
% Plot the mesh with curvature
trisurf(tri, pts(:,1), pts(:,2), pts(:,3), K, 'EdgeColor', 'none');
colormap jet;
colorbar;
title('Gaussian Curvature of Reconstructed Surface');
xlabel('X');
ylabel('Y');
zlabel('Z');
axis equal;
You can also refer to the following file exchange link for surface curvature computation which will help you with computations and plotting of Gaussian curvatures: https://www.mathworks.com/matlabcentral/fileexchange/11168-surface-curvature
