Merge pull request #4948 from nortikin/fix_4947_overlapping_voronoi_3…

…d_sites

#4947 fixing Overlapping Voronoi 3D sites, increase performance.

utils/geom.py

@@ -45,6 +45,29 @@
 from sverchok.dependencies import numba  # not strictly needed i think...
 from sverchok.utils.decorators_compilation import njit
 
+def bounding_box_aligned(verts):
+    # based on "3D Oriented bounding boxes": https://logicatcore.github.io/scratchpad/lidar/sensor-fusion/jupyter/2021/04/20/3D-Oriented-Bounding-Box.html
+    data = np.vstack(np.array(verts).transpose())
+    means = np.mean(data, axis=1)
+
+    cov = np.cov(data)
+    eval, evec = np.linalg.eig(cov)
+    centered_data = data - means[:,np.newaxis]
+    xmin, xmax, ymin, ymax, zmin, zmax = np.min(centered_data[0, :]), np.max(centered_data[0, :]), np.min(centered_data[1, :]), np.max(centered_data[1, :]), np.min(centered_data[2, :]), np.max(centered_data[2, :])
+    aligned_coords = np.matmul(evec.T, centered_data)
+    xmin, xmax, ymin, ymax, zmin, zmax = np.min(aligned_coords[0, :]), np.max(aligned_coords[0, :]), np.min(aligned_coords[1, :]), np.max(aligned_coords[1, :]), np.min(aligned_coords[2, :]), np.max(aligned_coords[2, :])
+
+    rectCoords = lambda x1, y1, z1, x2, y2, z2: np.array([[x1, x1, x2, x2, x1, x1, x2, x2],
+                                                        [y1, y2, y2, y1, y1, y2, y2, y1],
+                                                        [z1, z1, z1, z1, z2, z2, z2, z2]])
+
+    realigned_coords = np.matmul(evec, aligned_coords)
+    realigned_coords += means[:, np.newaxis]
+
+    rrc = np.matmul(evec, rectCoords(xmin, ymin, zmin, xmax, ymax, zmax))
+    rrc += means[:, np.newaxis]
+    rrc = rrc.transpose()
+    return tuple([rrc])
 
 identity_matrix = Matrix()
 

utils/voronoi3d.py

@@ -19,6 +19,7 @@
 import numpy as np
 from collections import defaultdict
 import itertools
+import datetime
 
 import bpy
 import bmesh
@@ -28,12 +29,12 @@
 from sverchok.data_structure import repeat_last_for_length
 from sverchok.utils.sv_mesh_utils import mask_vertices, polygons_to_edges, point_inside_mesh
 from sverchok.utils.sv_bmesh_utils import bmesh_from_pydata, pydata_from_bmesh, bmesh_clip
-from sverchok.utils.geom import calc_bounds, bounding_sphere, PlaneEquation
+from sverchok.utils.geom import calc_bounds, bounding_sphere, PlaneEquation, bounding_box_aligned
 from sverchok.utils.math import project_to_sphere, weighted_center
 from sverchok.dependencies import scipy, FreeCAD
 
 if scipy is not None:
-    from scipy.spatial import Voronoi, SphericalVoronoi
+    from scipy.spatial import Voronoi, SphericalVoronoi, Delaunay
 
 if FreeCAD is not None:
     from FreeCAD import Base
@@ -234,83 +235,179 @@ def calc_bvh_projections(bvh, sites):
             projections.append(loc)
     return np.array(projections)
 
-def voronoi_on_mesh_bmesh(verts, faces, n_orig_sites, sites, spacing=0.0, fill=True, precision=1e-8):
+# see additional info https://github.com/nortikin/sverchok/pull/4948
+def voronoi_on_mesh_bmesh(verts, faces, n_orig_sites, sites, spacing=0.0, mode='VOLUME', precision=1e-8):
 
-    def get_ridges_per_site(voronoi):
+    def get_sites_delaunay_params(delaunay):
         result = defaultdict(list)
-        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]
+        ridges = []
+        sites_pair = dict()
+        for simplex in delaunay.simplices:
+            ridges += itertools.combinations(tuple( sorted( simplex ) ), 2)
+
+        ridges = list(set( ridges )) # remove duplicates of ridges
+        ridges.sort() # for nice view in debugger
+
+        for ridge_idx in range(len(ridges)):
+            site1_idx, site2_idx = tuple(ridges[ridge_idx])
+            site1 = delaunay.points[site1_idx]
+            site2 = delaunay.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)
-        return result
+            normal =  Vector(site1 - site2).normalized() # normal to site1
+            plane1 = PlaneEquation.from_normal_and_point( normal, middle)
+            plane2 = PlaneEquation.from_normal_and_point(-normal, middle)
+            result[site1_idx].append( (site2_idx, site1, site2, middle,  normal, plane1) )
+            result[site2_idx].append( (site1_idx, site2, site1, middle, -normal, plane2) )
 
-    def cut_cell(verts, faces, planes, site, spacing):
-        src_mesh = bmesh_from_pydata(verts, [], faces, normal_update=True)
-        n_cuts = 0
-        for plane in planes:
-            if len(src_mesh.verts) == 0:
-                break
-            geom_in = src_mesh.verts[:] + src_mesh.edges[:] + src_mesh.faces[:]
-
-            plane_co = plane.projection_of_point(site)
-            plane_no = plane.normal.normalized()
-            if plane.side_of_point(site) > 0:
-                plane_no = - plane_no
-
-            plane_co = plane_co - 0.5 * spacing * plane_no
-
-            current_verts = np.array([tuple(v.co) for v in src_mesh.verts])
-            signs = PlaneEquation.from_normal_and_point(plane_no, plane_co).side_of_points(current_verts)
-            #print(f"Plane co {plane_co}, no {plane_no}, signs {signs}")
-            if (signs <= 0).all():# or (signs <= 0).all():
-                continue
+        return result
 
-            res = bmesh.ops.bisect_plane(
-                    src_mesh, geom=geom_in, dist=precision,
-                    plane_co = plane_co,
-                    plane_no = plane_no,
-                    use_snap_center = False,
-                    clear_outer = True,
-                    clear_inner = False
-                )
-            n_cuts += 1
-
-            if fill:
-                surround = [e for e in res['geom_cut'] if isinstance(e, bmesh.types.BMEdge)]
-                if surround:
-                    fres = bmesh.ops.edgenet_prepare(src_mesh, edges=surround)
-                    if fres['edges']:
-                        bmesh.ops.edgeloop_fill(src_mesh, edges=fres['edges'])
-
-        if n_cuts == 0:
+    # some statistics:
+    num_bisect = 0 # general count of bisect for full cutting process
+    num_unpredicted_erased = 0 # if optimisation can not find a skip bisect case (with using bounding box) then counter incremented
+
+    def cut_cell(start_mesh, sites_delaunay_params, site_idx, spacing, center_of_mass, bbox_aligned):
+        nonlocal num_bisect, num_unpredicted_erased
+        site_params = sites_delaunay_params[site_idx]
+
+        if len(start_mesh.verts) > 0:
+            lst_ridges_to_bisect = []
+            arr_dist_site_middle = np.empty(0)
+
+            out_of_bbox = False
+            src_mesh = None
+
+            # Sorting for optiomal bisections and search what can be skipped:
+            for i, (site_pair_idx, site_vert, site_pair_vert, middle, plane_no, plane) in enumerate(site_params):
+                # Move bisect plane on size of half of spacing (normal point to the site_idx from site_pair_idx)
+                plane_co = middle + 0.5 * spacing * plane_no
+                # [1]. Test if bbox_aligned outside a site_pair plane?
+                signs_verts_bbox_aligned = PlaneEquation.from_normal_and_point( plane_no, plane_co ).side_of_points(bbox_aligned)
+                # if all vertexes of bbox_aligned out of plane with negation normal then object will be erased anyway.
+                # So one can skeep bisect operation
+                if (signs_verts_bbox_aligned <= 0).all():
+                    out_of_bbox = True
+                    break
+                # if all vertexes of bbox_aligned is on a positive side of a plane then bisect cannot produce any sections.
+                # So one can skip operation of bisection and stay object unchanged (do not add ringe to bisection list)
+                if (signs_verts_bbox_aligned > 0).all():
+                    pass
+                else:
+                    # [2]. calc middle planes for optimal bisects sequence (sort later)
+                    plane_spacing = PlaneEquation.from_normal_and_point(plane_no, plane_co)
+                    sign = plane_spacing.side_of_points(center_of_mass)
+                    dist = plane_spacing.distance_to_point(center_of_mass)
+
+                    lst_ridges_to_bisect.append( [dist*sign, site_pair_idx, site_vert, site_pair_vert, middle, plane_co, plane_no, plane, ] )
+
+                # [3]. for test if all (site, middle) dist are less 0.5 spacing?
+                #    if spacing to big and eat all area [all (site-middle).lenght <= spacing/2]
+                arr_dist_site_middle = np.append(arr_dist_site_middle, np.linalg.norm(site_vert-middle) )
+
+                # here is the place to extend optimization variants to exclude bisect from process.
+                # To the future: one cannot optimize process of bisection. Only count of bisects can be optimized.
+                pass
+
+            # (3).
+            # out_of_bbox may realized before all site pairs observed so arr_dist_site_middle may contain not all dists
+            if out_of_bbox==False and (arr_dist_site_middle<=0.5 * spacing).all():
+                out_of_bbox = True
+                pass
+
+            if out_of_bbox==False:
+                # (2)
+                lst_ridges_to_bisect.sort()  # less dist gets more points to cut off (with negative dists to. Negative dist is a negative side of bisect plane)
+
+                src_mesh = start_mesh.copy() # do not need create src_mesh until here.
+
+                # A main bisection process of site_idx
+                for i in range(len(lst_ridges_to_bisect)):
+                    dist_center_of_mass_to_plane, site_pair_idx, site_vert, site_pair_vert, middle, plane_co, plane_no, plane = lst_ridges_to_bisect[i]
+                    geom_in = src_mesh.verts[:] + src_mesh.edges[:] + src_mesh.faces[:]
+                    res_bisect = bmesh.ops.bisect_plane(
+                            src_mesh, geom=geom_in, dist=precision,
+                            plane_co = plane_co,
+                            plane_no = plane_no,
+                            use_snap_center = False,
+                            clear_outer = False,
+                            clear_inner = True
+                        )
+                    num_bisect+=1 # for statistics
+
+                    if len(res_bisect['geom_cut'])>0:
+                        if mode=='VOLUME': # fill faces after bisect
+                            surround = [e for e in res_bisect['geom_cut'] if isinstance(e, bmesh.types.BMEdge)]
+                            if surround:
+                                fres = bmesh.ops.edgenet_prepare(src_mesh, edges=surround)
+                                if fres['edges']:
+                                    #bmesh.ops.edgeloop_fill(src_mesh, edges=fres['edges']) # has glitches
+                                    mfilled = bmesh.ops.triangle_fill(src_mesh, use_beauty=True, use_dissolve=True, edges=fres['edges'])
+                                else:
+                                    pass
+                            else:
+                                pass
+                    else:
+                        # if no geometry after bisect then break
+                        # Geometry get clear in two cases:
+                        # 1. Optimisation fail and not realized that this process has no result
+                        # 2. Big spacing eat geometry inside mesh
+                        if len( res_bisect['geom'] )==0:
+                            num_unpredicted_erased+=1 # for statistics
+                            break
+                        pass
+            else:
+                # func come here if out_of_bbox==True
+                pass
+
+        # if out_of_bbox==True then bisect process jump here
+        # if no verts then return noting
+        if src_mesh is None or len( src_mesh.verts ) == 0:
+            if src_mesh is not None:
+                src_mesh.clear() #remember to clear empty geometry
+                src_mesh.free()
             return None
 
-        return pydata_from_bmesh(src_mesh)
+        # if src_mesh has vertices then return mesh data
+        pydata = pydata_from_bmesh(src_mesh)
+        return pydata
 
     verts_out = []
     edges_out = []
     faces_out = []
 
-    voronoi = Voronoi(np.array(sites))
-    ridges_per_site = get_ridges_per_site(voronoi)
+    # http://www.qhull.org/html/qdelaun.htm
+    # http://www.qhull.org/html/qh-optc.htm
+    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html
+    delaunay = Delaunay(np.array(sites, dtype=np.float32))
+    sites_delaunay_params = get_sites_delaunay_params(delaunay)
+
     if isinstance(spacing, list):
         spacing = repeat_last_for_length(spacing, len(sites))
     else:
         spacing = [spacing for i in range(len(sites))]
+
+    # calc center of mass. Using for sort of bisect planes for sites.
+    center_of_mass = np.average( verts, axis=0 )
+    # using for precalc unneeded bisects
+    bbox_aligned = bounding_box_aligned(verts)[0]
+
+    start_mesh = bmesh_from_pydata(verts, [], faces, normal_update=False)
     for site_idx in range(len(sites)):
-        cell = cut_cell(verts, faces, ridges_per_site[site_idx], sites[site_idx], spacing[site_idx])
+        cell = cut_cell(start_mesh, sites_delaunay_params, site_idx, spacing[site_idx], center_of_mass, bbox_aligned)
         if cell is not None:
             new_verts, new_edges, new_faces = cell
             if new_verts:
                 verts_out.append(new_verts)
                 edges_out.append(new_edges)
                 faces_out.append(new_faces)
-
+    start_mesh.clear() # remember to clear empty geometry
+    start_mesh.free()
+
+
+    # show statistics:
+    # bisects - count of bisects in cut_cell
+    # unb - unpredicted erased mesh (bbox_aligned cannot make predicted results)
+    # sites - count of sites in process
+    # print( f"bisects: {num_bisect: 4d}, unb={num_unpredicted_erased: 4d}, sites={len(sites)}")
     return verts_out, edges_out, faces_out
 
 def voronoi_on_mesh(verts, faces, sites, thickness,
@@ -345,15 +442,8 @@ def voronoi_on_mesh(verts, faces, sites, thickness,
 
     else: # VOLUME, SURFACE
         all_points = sites[:]
-        if do_clip:
-            clipping = float(clipping)
-            for site in sites:
-                loc, normal, index, distance = bvh.find_nearest(site)
-                if loc is not None:
-                    p1 = loc + clipping * normal
-                    all_points.append(p1)
         verts, edges, faces = voronoi_on_mesh_bmesh(verts, faces, len(sites), all_points,
-                spacing = spacing, fill = (mode == 'VOLUME'),
+                spacing = spacing, mode = mode,
                 precision = precision)
         return verts, edges, faces, all_points