simple_delaunay_refinement_3d_for_blender.py

import bpy
import numpy as np
from scipy.spatial import ConvexHull

MAX_ITERATION = 10
MAX_ERROR_SQUARED = 1e-4
NUM_POINTS = 400
RADIUS = 10  # 球の半径

points = np.random.randn(3, NUM_POINTS)
points = (points / np.linalg.norm(points, axis=0)).T

faces = ConvexHull(points).simplices

for n in range(MAX_ITERATION):
    new_points = []  # 調整後の点のリスト
    ok = True
    max_error = -1

    # 各点(i=index)について以下を行う
    for i in range(len(points)):
        # 点i
        point = points[i]
        # 点iを含む三角形を抽出
        idxs = faces[(faces == i).any(axis=1)]
        # indexの重複を取り除き、各indexの出現回数を取得
        idxs, counts = np.unique(idxs, return_counts=1)
        # 点iは以下の計算から除外(...しなくてもよい？)
        idxs = idxs[idxs != i]
        # indexに対応する座標を取得
        ps = points.take(idxs, axis=0)
        # 座標の平均値を取得。調整後の点iの位置とする
        avg = np.mean(ps, axis=0)
        #if np.isnan(avg).any():
        #    continue
        avg = avg / np.linalg.norm(avg) # 正規化
        
        # 元の座標との距離(の２乗)が大きければ要再処理とする
        error = np.sum((avg - point)**2)
        if error > MAX_ERROR_SQUARED:
            max_error = max(error, max_error)
            ok = False
        new_points.append(avg)
        
    print(n, max_error)
        
    points = np.array(new_points) 
    faces = ConvexHull(points).simplices
    if ok:
        break
    
print("create a mesh")
points *= RADIUS
mesh = bpy.data.meshes.new("sphere")
obj = bpy.data.objects.new(mesh.name, mesh)
col = bpy.data.collections.get("Collection")
col.objects.link(obj)
bpy.context.view_layer.objects.active = obj

mesh.from_pydata(points.tolist(), [], faces.tolist())