Curvature computation

[frk1993]

Hi guys, Is there a function/method to compute curvature of the points?

Thanks

[griegler]

No, at the moment there is no method for this. Would you be interested in contributing one with our help?

[frk1993]

Thank you for the quick response.

Yes, it would be great.

Thanks!

[nickponline]

I would also be very interested in this.

[nickponline]

@frk1993 @griegler here is a way to do it although slow, any comments of the preferred way open3d would handle this?

def compute_curvature(pcd, radius=0.5):

    points = np.asarray(pcd.points)

    from scipy.spatial import KDTree
    tree = KDTree(points)

    curvature = [ 0 ] * points.shape[0]

    for index, point in enumerate(points):
        indices = tree.query_ball_point(point, radius)

        # local covariance
        M = np.array([ points[i] for i in indices ]).T
        M = np.cov(M)

        # eigen decomposition
        V, E = np.linalg.eig(M)
        # h3 < h2 < h1
        h1, h2, h3 = V

        curvature[index] = h3 / (h1 + h2 + h3)

    return curvature

[griegler]

It should be very simple to add such a function to Open3D, should be pretty similar to the EstimatesNormals here.

[MehdiMaboudi]

Any updates on this topic?

[ljmartin]

I'm happy to try this, but the function above is not giving results that look like curvature to me. is that how it should look? Areas with opposite curvature seem to be the same colour. Does anyone know a reference paper or implementation?

import open3d as o3d
import meshplot as mp
armadillo_mesh = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo_mesh.path)

curvatures = compute_curvature(np.array(mesh.vertices), radius=5)

mp.plot(np.array(mesh.vertices),
        np.array(mesh.triangles), c=np.array(curvatures))

:

[tillns]

I'm happy to try this, but the function above is not giving results that look like curvature to me. is that how it should look? Areas with opposite curvature seem to be the same colour. Does anyone know a reference paper or implementation?

import open3d as o3d
import meshplot as mp
armadillo_mesh = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo_mesh.path)

curvatures = compute_curvature(np.array(mesh.vertices), radius=5)

mp.plot(np.array(mesh.vertices),
        np.array(mesh.triangles), c=np.array(curvatures))

:

I'm currently using IGL's implementation for computing the principal curvatures. This is their implementation: code link with the doc here: doc link. Would be cool to have something similar in open3d.

[AlexWang1900]

How about this:

import open3d as o3d
import numpy as np

def pca_compute(data, sort=True):
    """
	SVD decomposition
    """
    average_data = np.mean(data, axis=0) 
    decentration_matrix = data - average_data
    H = np.dot(decentration_matrix.T, decentration_matrix)
    eigenvectors, eigenvalues, eigenvectors_T = np.linalg.svd(H)
    if sort:
        sort = eigenvalues.argsort()[::-1]
        eigenvalues = eigenvalues[sort]
    return eigenvalues

def caculate_surface_curvature(cloud, radius=0.003):
    points = np.asarray(cloud.points)
    kdtree = o3d.geometry.KDTreeFlann(cloud)
    num_points = len(cloud.points)
    curvature = []  
    for i in range(num_points):
        k, idx, _ = kdtree.search_radius_vector_3d(cloud.points[i], radius)
        neighbors = points[idx, :]
        w = pca_compute(neighbors)
        delt = np.divide(w[2], np.sum(w), out=np.zeros_like(w[2]), where=np.sum(w) != 0)
        curvature.append(delt)
    curvature = np.array(curvature, dtype=np.float64)
    return curvature

if __name__ == ''__main__'':
    pcd = o3d.io.read_point_cloud("data/1.pcd")
    surface_curvature = caculate_surface_curvature(pcd, radius=0.003)
    print(surface_curvature[:10]) 

[chibai]

No offense, but any update in this topic??

[manhha1402]

Will be there the update for curvature of pointcloud?

[parsazarei1]

@frk1993 @griegler here is a way to do it although slow, any comments of the preferred way open3d would handle this?

def compute_curvature(pcd, radius=0.5):

    points = np.asarray(pcd.points)

    from scipy.spatial import KDTree
    tree = KDTree(points)

    curvature = [ 0 ] * points.shape[0]

    for index, point in enumerate(points):
        indices = tree.query_ball_point(point, radius)

        # local covariance
        M = np.array([ points[i] for i in indices ]).T
        M = np.cov(M)

        # eigen decomposition
        V, E = np.linalg.eig(M)
        # h3 < h2 < h1
        h1, h2, h3 = V

        curvature[index] = h3 / (h1 + h2 + h3)

    return curvature

IS this gaussian or mean curvature? how would you calculate both ?

[parsazarei1]

@nickponline

[AlbertoMQ]

Trimesh has an implementation for different kinds of curvature. https://github.com/mikedh/trimesh

[parsazarei1]

Hi Alberto, Thanks for sending me this! I use trimesh to find the curvature.

Best,

From: Alberto Miller @.> Sent: Thursday, February 23, 2023 4:04 PM To: isl-org/Open3D @.> Cc: Parsa Zareiesfandabadi @.>; Comment @.> Subject: Re: [isl-org/Open3D] Curvature computation (#2471)

Trimesh has an implementation for different kinds of curvature. https://github.com/mikedh/trimeshhttps://urldefense.com/v3/__https:/github.com/mikedh/trimesh__;!!OToaGQ!vBAnxGgAhr18mo2tmZtMu_Lceh4zFBDd-t7QnwnYJwYHyNtoCzFi6M7mXemiUyK2fWbcZ75UkZYCZ7IIjEabuKik$

— Reply to this email directly, view it on GitHubhttps://urldefense.com/v3/__https:/github.com/isl-org/Open3D/issues/2471*issuecomment-1442428463__;Iw!!OToaGQ!vBAnxGgAhr18mo2tmZtMu_Lceh4zFBDd-t7QnwnYJwYHyNtoCzFi6M7mXemiUyK2fWbcZ75UkZYCZ7IIjDNSFnYd$, or unsubscribehttps://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/A5U6FSKJRRHYIPOKXVDWWCTWY7GC3ANCNFSM4SQWF5YQ__;!!OToaGQ!vBAnxGgAhr18mo2tmZtMu_Lceh4zFBDd-t7QnwnYJwYHyNtoCzFi6M7mXemiUyK2fWbcZ75UkZYCZ7IIjLEPKag7$. You are receiving this because you commented.Message ID: @.@.>>

[ssheorey]

Open3D now has covariance computation.

http://www.open3d.org/docs/latest/python_api/open3d.geometry.PointCloud.html#open3d.geometry.PointCloud.estimate_covariances

This looks like a good first issue if someone wants to contribute curvature computation.

[HeCraneChen]

I recently developed a new tool for total curvature estimation. The code compiles and runs on MacOS, Ubuntu, and Windows. Feel free to grab the tool here:

https://github.com/HeCraneChen/total-curvature-estimation.git

[HeCraneChen]

The tool I developed has a open3D version as well (compilation tested just on MacOS so far). Feel free to grab the open3d-style code here: https://github.com/HeCraneChen/open3d-discrete-total-curvature.git

[HeCraneChen]

Open3D now has covariance computation.

http://www.open3d.org/docs/latest/python_api/open3d.geometry.PointCloud.html#open3d.geometry.PointCloud.estimate_covariances

This looks like a good first issue if someone wants to contribute curvature computation.

I'm interested in contributing curvature computation, let me know if you have any leads ;) Thinkings about covariance matrix: It can be reasonable if the point cloud is uniformly sampled. However, the normalization process might pose a problem if the distribution of the point cloud is irregular.

[yyvhang]

The tool I developed has a open3D version as well (compilation tested just on MacOS so far). Feel free to grab the open3d-style code here: https://github.com/HeCraneChen/open3d-discrete-total-curvature.git

Does this have a python version?

[zxhcho]

I'm happy to try this, but the function above is not giving results that look like curvature to me. is that how it should look? Areas with opposite curvature seem to be the same colour. Does anyone know a reference paper or implementation?

import open3d as o3d
import meshplot as mp
armadillo_mesh = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo_mesh.path)

curvatures = compute_curvature(np.array(mesh.vertices), radius=5)

mp.plot(np.array(mesh.vertices),
        np.array(mesh.triangles), c=np.array(curvatures))

:

I'm currently using IGL's implementation for computing the principal curvatures. This is their implementation: code link with the doc here: doc link. Would be cool to have something similar in open3d.

The returned eigenvalue of np.linalg.eig() is not supposed to be ordered, so this code should be curvature[index] = min(V)/sum(V)

[HeCraneChen]

The tool I developed has a open3D version as well (compilation tested just on MacOS so far). Feel free to grab the open3d-style code here: https://github.com/HeCraneChen/open3d-discrete-total-curvature.git

Does this have a python version?

Not yet. I indeed plan to work on a python wrapper, but have been busy, therefore it is not ready yet. I assure you though, the C++ version is clean and easy to setup/compile/run. If you want to give the C++ version a shot, I'm happy to answer all your questions if you ran into any compilation problems.

[pnchuyen]

Here is faster based on the above suggestions:

def calculate_surface_curvature(pcd, radius=0.1, max_nn=30):
    pcd_n = copy.deepcopy(pcd)
    pcd_n.estimate_covariances(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=radius, max_nn=max_nn))
    covs = np.asarray(pcd_n.covariances)
    vals, vecs = np.linalg.eig(covs)
    curvature = np.min(vals, axis=1)/np.sum(vals, axis=1)
    return curvature