How can I calculate the gussian curvature of a point cloud by matlab?

Hello every body. Does exist any matlab code or built in function for gussian curvature computatation by matlab?

Answer 1

[Divyam]

Hi @Ashkan Rigi,

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