could we support python usage？

[Fil]

I don't plan to do that, but you're absolutely welcome to port this to python. The underlying library (https://github.com/mapbox/delaunator) already has a python port. Please comment here if you have questions about the algo (and to report success!).

[MengHao666]

https://github.com/mapbox/delaunator Thanks for your warm reply! Your work is very impressive!

Aactually, I want to get convex or non-convex hull from points on surface of the unit sphere. the recommended repo (https://github.com/mapbox/delaunator) could only support 2d points.

The feature I want is [voronoi.hull(data)](https://github.com/Fil/d3-geo-voronoi/blob/main/src/voronoi.js）. It is a really exciting feature for me . and I read some source code try to re-implement. I try to do Sphere Delaunay triangulation first and then compute the edges which is used for one triangle only to define the boundary hull of these points.

However, I failed. I just use scipy.spatial.Delaunay and [SphericalVoronoi]https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.SphericalVoronoi.html#scipy.spatial.SphericalVoronoi as i found the desciption is quite similar. It might be that I did not know the related math knowledge. And I felt quite upset for that.

Finally, I found my data is quite special. it is the rotation of a finger joint, so it is restricted in a hemisphere. So I just find a plane and project ponits into that, and get 2d convex hull indices, and finally remap it into sphere. It is quite dirty but seems works.

However, I will need to compute if the point in such sphere hull and the minimal distance from a point to the boundary of the sphere hull.

I will be very thankful if you could give some mathmatical hints for my wrong usgae of scipy.spatial.Delaunay or SphericalVoronoi. Perhaps, you could give me some further explaination about your source code, i tried to translate it into python with gpt4 but it looks failed.

thanks, hao

[Fil]

To compute the hull on a hemisphere, you certainly don't need all the complexities involved in spherical delaunay. Transform that hemisphere to the plane with the stereographic projection, and take the hull from the projected points. The crucial mathematical property (for this task) of the stereographic is that it maps circles to circles, and thus preserves the Delaunay property.

[MengHao666]

To compute the hull on a hemisphere, you certainly don't need all the complexities involved in spherical delaunay. Transform that hemisphere to the plane with the stereographic projection, and take the hull from the projected points. The crucial mathematical property (for this task) of the stereographic is that it maps circles to circles, and thus preserves the Delaunay property.

Thanks for your warm reply. In my own situation, the point on unit sphere is restircted in hemisphere. I mean it might not restrict in the north part or south part, but in the western part or east one.

I search a lot in the internet about "stereographic projection,". And I found it might be mainly realted to geography plotting map. I am struggled in how to use or understand it.

my current result is like that

my current result is computed from following code with the very simple projection, "vectors" is just the point on surface:

import math
import os

import numpy as np
import alphashape
from matplotlib import pyplot as plt
from scipy.spatial import ConvexHull
from scipy.spatial import geometric_slerp
from tqdm import tqdm

 from mobu.scripts.python39.MocapRefinement.modules.cleanup_finger.utils.rot_util import \
      cartesian_to_spherical, \
      get_init_bone_vector_rot_meaning, swing_twist_decomposition, angle_between_vectors, find_beta_gamma, \
      vectors2spherical, vectors2hull

  def projection_to_plane(point, normal):
      return point - np.dot(point, normal) * normal

  def projection(point, direction):
      return np.dot(point, direction)

  def get_concave_hull(points_2d, alpha_value=0):
      point_to_index = {tuple(point): idx for idx, point in enumerate(points_2d)}

      alpha_shape = alphashape.alphashape(points_2d, alpha_value)
      x, y = alpha_shape.exterior.xy

      ordered_points = np.zeros([len(x), 2])
      ordered_points[:, 0] = x
      ordered_points[:, 1] = y

      indices = [point_to_index[point] for point in alpha_shape.exterior.coords]

      return ordered_points, indices

  def find_plane(points):
      if len(points) < 3:
          return "需要至少3个点来确定一个平面"

      max_distance_sum = -np.inf
      best_direction = None

      theta_values = np.linspace(0, 2 * np.pi, 100)
      phi_values = np.linspace(0, np.pi, 100)

      theta_matrix, phi_matrix = np.meshgrid(theta_values, phi_values)
      direction_matrix = np.array([np.sin(phi_matrix) * np.cos(theta_matrix),
                                   np.sin(phi_matrix) * np.sin(theta_matrix),
                                   np.cos(phi_matrix)]).T

      for direction in direction_matrix.reshape(-1, 3):
          projections = np.dot(points, direction)
          distance_sum = np.sum(np.abs(projections))

          if distance_sum > max_distance_sum and (np.all(projections >= 0) or np.all(projections <= 0)):
              max_distance_sum = distance_sum
              best_direction = direction

      return best_direction

   def visualize(vectors, normal, projections, hull_indices, title):
    # plot 平面，法向量
    fig = plt.figure(figsize=(20, 15))

    # 绘制单位球
    radius = 1
    # ax1 = fig.add_subplot(121, projection='3d')
    ax1 = fig.add_subplot(projection='3d')
    u = np.linspace(0, 2 * np.pi, 100)  # 用参数方程画图
    v = np.linspace(0, np.pi, 100)
    x = radius * np.outer(np.cos(u), np.sin(v))  # outer()外积函数：返回cosu × sinv
    y = radius * np.outer(np.sin(u), np.sin(v))  #
    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))  # ones()：返回一组[1,1,.......]
    ax1.plot_surface(x, y, z, color='g', alpha=0.3)

    # 绘制所有的bone的点
    ax1.scatter(*(vectors).T, color='r', s=1)

    # 绘制投影点
    # ax1.scatter(*(projections).T, color='b', s=1)

    #绘制求得的hull
    # for s in hull_indices:
    #     ax1.plot(projections[s, 0], projections[s, 1], projections[s, 2], "k-")

    # plot 平面
    # if isinstance(normal, np.ndarray):
    #     x_range = np.linspace(-1, 1, 10)
    #     y_range = np.linspace(-1, 1, 10)
    #     X, Y = np.meshgrid(x_range, y_range)
    #     Z = (-normal[0] * X - normal[1] * Y) / normal[2]
    #     ax1.plot_surface(X, Y, Z, alpha=1)
    # ax1.scatter(*(stance_bone_vector).T, marker='P', color='y', s=200)
    # ax1.scatter(1, 1,
    #            1, color='r', s=100)

    # hull
    # print(hull.simplices)
    # sphere_hull = vectors[hull.simplices]
    t_vals = np.linspace(0, 1, 2000)
    # for i,j in convex_hull_indices:
    for i, j in hull_indices:
        start = vectors[i]
        end = vectors[j]
        # print(start)
        # print(end)
        result = geometric_slerp(start, end, t_vals)
        ax1.plot(result[..., 0],
                 result[..., 1],
                 result[..., 2],
                 c='y',
                 linewidth=2, markersize=5)

    ax1.set_xlabel("x'")
    ax1.set_ylabel("y'")
    ax1.set_zlabel("z'")
    ax1.set_aspect('equal')
    ax1.set_title("旋转空间")

    plt.grid(True)
    plt.grid(color='r', linestyle='--', linewidth=1, alpha=0.3)

    # ax2 = fig.add_subplot(122)
    # ax2.scatter(points_2d[:, 0], points_2d[:, 1], s=20)
    # ax2.set_xlabel("x'")
    # ax2.set_ylabel("y'")
    # plt.grid(True)
    # plt.grid(color='r', linestyle='--', linewidth=1, alpha=0.3)
    # for s in hull.simplices:
    #     ax2.plot(points_2d[s, 0], points_2d[s, 1], "k-")

    plt.suptitle(title)

    plt.show()

   def process_vectors(vectors, txt_name):
      normal = find_plane(vectors)
      # print("法向量:", normal)

      if not isinstance(normal, np.ndarray):
          raise Exception(f"")

      # # 绘制平面
      # if isinstance(normal, np.ndarray):
      # x_range = np.linspace(-1, 1, 10)
      # y_range = np.linspace(-1, 1, 10)
      # X, Y = np.meshgrid(x_range, y_range)
      # Z = (-normal[0] * X - normal[1] * Y) / normal[2]
      # ax1.plot_surface(X, Y, Z, alpha=1)

      # 计算投影点并绘制
      projections = np.array([projection_to_plane(point, normal) for point in vectors])
      # print(projections.shape)

      # 将投影点进行转化
      # 选择一个参考点作为原点
      new_origin = projections[0]
      # 计算x'轴和y'轴
      new_x_axis = np.array([1, 0, 0]) if np.dot(normal, np.array([1, 0, 0])) < 0.9 else np.array([0, 1, 0])
      new_y_axis = np.cross(normal, new_x_axis)

      # 计算每个点在新坐标系中的坐标
      new_points_2d = np.array(
          [[(p - new_origin).dot(new_x_axis), (p - new_origin).dot(new_y_axis)] for p in projections])

      # 计算凸包并绘制
      convex_hull = ConvexHull(new_points_2d, qhull_options="QJ")
      convex_hull_indices = convex_hull.simplices

      # 计算凹包
      concave_hull_ordered_points, concave_hull_order_indices = get_concave_hull(new_points_2d, alpha_value=2.0)
      concave_hull_indices = [[concave_hull_order_indices[i], concave_hull_order_indices[i + 1]] for i in
                              range(len(concave_hull_order_indices)) if i < len(concave_hull_order_indices) - 1]

      visualize(vectors, normal, projections, convex_hull_indices, title=f"{txt_name} - convex_hull")
      visualize(vectors, normal, projections, concave_hull_indices, title=f"{txt_name} - concave_hull")

[Fil]

I'm sorry I can't help with this.

[MengHao666]

I'm sorry I can't help with this.

thank you the same. Thank u again for your excellent work!