voronoi3d.py

"""
in sites_in v
in verts_in v
in faces_in s
out verts_out v
out faces_out s
"""

import numpy as np
from collections import defaultdict

import bmesh

from sverchok.data_structure import zip_long_repeat
from sverchok.utils.geom import linear_approximation, PlaneEquation
from sverchok.utils.sv_bmesh_utils import pydata_from_bmesh, bmesh_from_pydata
from sverchok.dependencies import scipy

if scipy is None:
    raise Exception("This script requires SciPy")

from scipy.spatial import Voronoi

def get_ridges_per_site(voronoi):
    result = defaultdict(list)
    print(voronoi.ridge_points)
    for ridge_idx in range(len(voronoi.ridge_points)):
        site1_idx, site2_idx = tuple(voronoi.ridge_points[ridge_idx])
        site1 = voronoi.points[site1_idx]
        site2 = voronoi.points[site2_idx]
        middle = (site1 + site2) * 0.5
        normal = site2 - site1
        plane = PlaneEquation.from_normal_and_point(normal, middle)
        result[site1_idx].append(plane)
        result[site2_idx].append(plane)
        if site1_idx in {27,22} and site2_idx in {27,22}:
            print(f"Ridge: #{site1_idx} - #{site2_idx} => {plane}")
    return result

def cut_cell(verts, faces, planes, site):
    src_mesh = bmesh_from_pydata(verts, [], faces, normal_update=True)
    for plane in planes:
        geom_in = src_mesh.verts[:] + src_mesh.edges[:] + src_mesh.faces[:]
        
        plane_co = plane.projection_of_point(site)
        plane_no = plane.normal
        if plane.side_of_point(site) > 0:
            plane_no = - plane_no
        #print(f"Plane co {plane_co}, no {plane_no}")
        res = bmesh.ops.bisect_plane(
                src_mesh, geom=geom_in, dist=0.00001,
                plane_co = plane_co,
                plane_no = plane_no,
                use_snap_center = False,
                clear_outer = True,
                clear_inner = False
            )
            
        surround = [e for e in res['geom_cut'] if isinstance(e, bmesh.types.BMEdge)]
        fres = bmesh.ops.edgenet_prepare(src_mesh, edges=surround)
        bmesh.ops.edgeloop_fill(src_mesh, edges=fres['edges'])  
    
    return pydata_from_bmesh(src_mesh)

verts_out = []
faces_out = []
for sites, verts, faces in zip_long_repeat(sites_in, verts_in, faces_in):
    #src_mesh = bmesh_from_pydata(verts, [], faces, normal_update=True)
    voronoi = Voronoi(np.array(sites))
    ridges_per_site = get_ridges_per_site(voronoi)
    for site_idx in range(len(sites)):
        new_verts, _, new_faces = cut_cell(verts, faces, ridges_per_site[site_idx], sites[site_idx])
        verts_out.append(new_verts)
        faces_out.append(new_faces)