From 098e32bf814d205a40acd4119d6bfc6540aed2c9 Mon Sep 17 00:00:00 2001 From: satabol Date: Tue, 3 Oct 2023 22:17:47 +0300 Subject: [PATCH] fix #5010. Refactor delaunay_3d_mk2 to automatically select topology - volumes, planes or edges. Full backward compatibility. --- nodes/spatial/delaunay_3d_mk2.py | 299 ++++++++++++++++++++++++++++--- 1 file changed, 274 insertions(+), 25 deletions(-) diff --git a/nodes/spatial/delaunay_3d_mk2.py b/nodes/spatial/delaunay_3d_mk2.py index e3de0a6f40..f131f53522 100644 --- a/nodes/spatial/delaunay_3d_mk2.py +++ b/nodes/spatial/delaunay_3d_mk2.py @@ -6,6 +6,7 @@ # License-Filename: LICENSE import numpy as np +from collections import defaultdict from itertools import combinations import bpy @@ -13,11 +14,145 @@ from sverchok.node_tree import SverchCustomTreeNode from sverchok.data_structure import updateNode, zip_long_repeat, ensure_nesting_level, get_data_nesting_level +from sverchok.utils.geom import PlaneEquation from sverchok.dependencies import scipy +from mathutils import Vector +import random if scipy is not None: from scipy.spatial import Delaunay + from scipy.spatial.transform import Rotation as R +def is_volume_0(verts, idxs, threshold): + if threshold == 0: + return False + a, b, c, d = [verts[i] for i in idxs] + a, b, c, d = np.array(a), np.array(b), np.array(c), np.array(d) + v1 = b - a + v2 = c - a + v3 = d - a + v1 = v1 / np.linalg.norm(v1) + v2 = v2 / np.linalg.norm(v2) + v3 = v3 / np.linalg.norm(v3) + volume = np.cross(v1, v2).dot(v3) / 6 + return abs(volume) < threshold + + +def is_area_0(verts, idxs, threshold): + if threshold == 0: + return False + a, b, c = [verts[i] for i in idxs] + a, b, c = np.array(a), np.array(b), np.array(c) + v1 = b - a + v2 = c - a + v1 = v1 / np.linalg.norm(v1) + v2 = v2 / np.linalg.norm(v2) + area = np.linalg.norm(np.cross(v1, v2)) / 2 + return abs(area) < threshold + +def is_length_0(verts, idxs, threshold): + if threshold == 0: + return False + a, b = [verts[i] for i in idxs] + a, b = np.array(a), np.array(b) + v1 = b - a + len = np.linalg.norm(v1) + return abs(len) < threshold + +def get_sites_delaunay_params(np_sites, n_orig_sites, threshold=0): + + delaunay = Delaunay(np_sites) + dimension_size_1 = -1 # simplex size on start + dimension_size_2 = -1 # simplex size of finish + result = defaultdict(list) + ridges = [] + simplices = [] + dict_source_dsimplex = dict() + sites_pair = dict() + for i, simplex in enumerate(delaunay.simplices): + dsimplex = tuple(simplex[simplex0: + # test simplex size for threshold. Remove simplexes that do not allowed in their dimensions. (plane is disallowed in a volume, line in disallowed in a plane). + arr = [] + vertices = np.array([delaunay.points[i,:3] for i,e in enumerate(delaunay.points)]) # point come with 4d vector coordinates + for sim in simplices: + simplex_size = len(sim) + # simplex is line. If line has size(length)==0 then do not use this simplex + if simplex_size==2 and is_length_0(vertices, sim, threshold): + continue + # simplex is plane. If plane has size(area)==0 then do not use this simplex + if simplex_size==3 and is_area_0(vertices, sim, threshold): + continue + # simplex is volumetric. If volume has size(volume)==0 then do not use this simplex + if simplex_size==4 and is_volume_0(vertices, sim, threshold): + continue + if simplex_size>4: + print(f"Very strange situation. size of simplex({simplex})=={simplex_size}. Skipped. Need attension") + continue + arr.append(sim) + pass + else: + arr = simplices + + # get length of first elem and select only elems with that len. + # this is a filter of allowed dimension of simplices. Top elem is a first elem of allowed dimensions. + arr_len = [] + len_0 = len(arr[0]) + dimension_size_2 = len_0 + for sim in arr: + if len(sim)==len_0: + arr_len.append(sim) + else: + break + list_removed_inner_simplices = arr_len + + #res = {"simplices": res} + list_restored_orders = [] + for key in list_removed_inner_simplices: + list_restored_orders.append( dict_source_dsimplex[tuple(key)] ) + #res = np.array(res.tolist(), dtype=np.int32) + res = np.array(list_restored_orders, dtype=np.int32) + return res, dimension_size_1, dimension_size_2 + + for ridge_idx in range(len(ridges)): + site1_idx, site2_idx = tuple(ridges[ridge_idx]) + # Remove 4D simplex ridges: + if n_orig_sites<=site1_idx or n_orig_sites<=site2_idx: + continue + # Convert source sites to the 3D + site1 = delaunay.points[site1_idx] + site1 = Vector([site1[0], site1[1], site1[2], ]) + site2 = delaunay.points[site2_idx] + site2 = Vector([site2[0], site2[1], site2[2], ]) + middle = (site1 + site2) * 0.5 + 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) ) + + return result class SvDelaunay3dMk2Node(SverchCustomTreeNode, bpy.types.Node): """ @@ -72,20 +207,6 @@ def make_faces(self, idxs): def get_verts(self, verts, idxs): return [verts[i] for i in idxs] - def is_planar(self, verts, idxs, threshold): - if threshold == 0: - return False - a, b, c, d = [verts[i] for i in idxs] - a, b, c, d = np.array(a), np.array(b), np.array(c), np.array(d) - v1 = b - a - v2 = c - a - v3 = d - a - v1 = v1 / np.linalg.norm(v1) - v2 = v2 / np.linalg.norm(v2) - v3 = v3 / np.linalg.norm(v3) - volume = np.cross(v1, v2).dot(v3) / 6 - return abs(volume) < threshold - def is_too_long(self, verts, idxs, threshold): if threshold == 0: return False @@ -121,16 +242,110 @@ def process(self): faces_item = [] for vertices, volume_threshold, edge_threshold in zip_long_repeat(*params): - tri = Delaunay(np.array(vertices)) + # https://github.com/nortikin/sverchok/pull/4952 + # 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 + # Convert sites to 4D + #np_sites = np.array(vertices, dtype=np.float32) + + # Think about oriented BBox: + # oX_size = np.max(np.array(vertices, dtype=np.float16)[:,0]) - np.min(np.array(vertices, dtype=np.float16)[:,0]) + # oY_size = np.max(np.array(vertices, dtype=np.float16)[:,1]) - np.min(np.array(vertices, dtype=np.float16)[:,1]) + # oZ_size = np.max(np.array(vertices, dtype=np.float16)[:,2]) - np.min(np.array(vertices, dtype=np.float16)[:,2]) + # axis = np.argsort([oX_size, oY_size, oZ_size]) + # plane_axis_X = axis[0] + # plane_axis_Y = axis[1] + # plane_axis_Z = axis[2] + # np_sites = np.array([(v[plane_axis_X], v[plane_axis_Y], v[plane_axis_Z], 0) for v in vertices], dtype=np.float32) + + np_sites = np.array([(*v, 0.0 ) for v in vertices], dtype=np.float32) + #np_sites = np.array( vertices, dtype=np.float32) + #np_sites = np.pad( np_sites, (0,1), 'constant') + #np_sites4 = np.zeros((len(np_sites), 4)) + #np_sites4[:,:3] = np_sites + for i in range(10): # 2 is max, but for insurance + # Add 3D tetraedre to the 4D with W=1 + np_sites = np.append(np_sites, [[0.0, 0.0, 0.0, 1.0], + [1.0, 0.0, 0.0, 1.0], + [0.0, 1.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 1.0], + ], axis=0) + + #delaunay = Delaunay(np.array(np_sites, dtype=np.float32)) + simplices, dim1, dim2 = get_sites_delaunay_params( np_sites, len(vertices), threshold=self.volume_threshold ) + if dim1==dim2: # 4 - volume, 3 - planes, 2 - edges, 1 - dots + print(f"Found solution for dim={dim2-1} with {i} attempt") + break + else: + oX_size = np.max(np.array(vertices, dtype=np.float16)[:,0]) - np.min(np.array(vertices, dtype=np.float16)[:,0]) + oY_size = np.max(np.array(vertices, dtype=np.float16)[:,1]) - np.min(np.array(vertices, dtype=np.float16)[:,1]) + oZ_size = np.max(np.array(vertices, dtype=np.float16)[:,2]) - np.min(np.array(vertices, dtype=np.float16)[:,2]) + axis = list(np.argsort([oX_size, oY_size, oZ_size])) + axis.reverse() + plane_axis_X = axis[0] + plane_axis_Y = axis[1] + plane_axis_Z = axis[2] + if dim2==3: # plane. faces and edges + # Find max plane axises of Bounding Box get it to XY plane + # convert vertices to plane oXY. get two first elems. Is is plane max area + np_sites = np.array([(v[plane_axis_X], v[plane_axis_Y], 0, 0) for v in vertices], dtype=np.float32) + continue + elif dim2==2: # line, only edges + # Find max side axis of Bounding Box get it to oX plane + np_sites = np.array([(v[plane_axis_X], 0, 0, 0) for v in vertices], dtype=np.float32) + continue + elif dim2==2: # dot. no edges + break + else: + # unknown dim2. insurence + break + else: + #print("delaunay_3d_mk2. Solution not found") # incredibly + pass + + #tri = Delaunay(np.array(vertices)) + #simplices = tri.simplices if self.join: verts_new = vertices edges_new = set() faces_new = set() - for simplex_idx, simplex in enumerate(tri.simplices): - if self.is_too_long(vertices, simplex, edge_threshold) or self.is_planar(vertices, simplex, volume_threshold): + for simplex_idx, simplex in enumerate(simplices): + # unknown simplex. skip. Incredible. get_sites_delaunay_params work with size==5 so need test. + if simplex.size>4: continue - edges_new.update(set(self.make_edges(simplex))) - faces_new.update(set(self.make_faces(simplex))) + + if self.is_too_long(vertices, simplex, edge_threshold): + continue + + # simplex is line. If line has size(length)==0 then do not use this simplex + if simplex.size==2 and is_length_0(vertices, simplex, volume_threshold): + continue + # simplex is plane. If plane has size(area)==0 then do not use this simplex + if simplex.size==3 and is_area_0(vertices, simplex, volume_threshold): + continue + # simplex is volumetric. If volume has size(volume)==0 then do not use this simplex + if simplex.size==4 and is_volume_0(vertices, simplex, volume_threshold): + continue + + # with simplex.size>=2 (2,3,4) then create edges + if simplex.size>=2: + #edges_simplex = self.make_edges(list(range(simplex.size))) + edges_simplex = self.make_edges(simplex) + edges_new.update(edges_simplex) + pass + + # with simplex.size>=3 (3,4) then create faces + if simplex.size>=3: + #faces_simplex = self.make_faces(list(range(simplex.size))) + faces_simplex = self.make_faces(simplex) + faces_new.update(faces_simplex) + pass + + # if some geometry get visible then geometery need verts: + #verts_simplex = self.get_verts(vertices, simplex) + #verts_new.extend(verts_simplex) + verts_item.append(verts_new) edges_item.append(list(edges_new)) faces_item.append(list(faces_new)) @@ -138,15 +353,48 @@ def process(self): verts_new = [] edges_new = [] faces_new = [] - for simplex in tri.simplices: - if self.is_too_long(vertices, simplex, edge_threshold) or self.is_planar(vertices, simplex, volume_threshold): + for simplex in simplices: + # unknown simplex. skip. Incredible. get_sites_delaunay_params work with size==5 so need test. + if simplex.size>4: continue + + if self.is_too_long(vertices, simplex, edge_threshold): + continue + + # simplex is line. If line has size(length)==0 then do not use this simplex + if simplex.size==2 and is_length_0(vertices, simplex, volume_threshold): + continue + # simplex is plane. If plane has size(area)==0 then do not use this simplex + if simplex.size==3 and is_area_0(vertices, simplex, volume_threshold): + continue + # simplex is volumetric. If volume has size(volume)==0 then do not use this simplex + if simplex.size==4 and is_volume_0(vertices, simplex, volume_threshold): + continue + + # with simplex.size>=2 (2,3,4) then create edges + if simplex.size>=2: + edges_simplex = self.make_edges(list(range(simplex.size))) + edges_new.append(edges_simplex) + else: + edges_new.append([]) + pass + + # with simplex.size>=3 (3,4) then create faces + if simplex.size>=3: + faces_simplex = self.make_faces(list(range(simplex.size))) + faces_new.append(faces_simplex) + else: + faces_new.append([]) + pass + + #edges_simplex = self.make_edges([0, 1, 2, 3]) + #faces_simplex = self.make_faces([0, 1, 2, 3]) + + # if some geometry get visible then geometery need verts: verts_simplex = self.get_verts(vertices, simplex) - edges_simplex = self.make_edges([0, 1, 2, 3]) - faces_simplex = self.make_faces([0, 1, 2, 3]) verts_new.append(verts_simplex) - edges_new.append(edges_simplex) - faces_new.append(faces_simplex) + + verts_item.extend(verts_new) edges_item.extend(edges_new) faces_item.extend(faces_new) @@ -158,6 +406,7 @@ def process(self): else: verts_out.extend(verts_item) edges_out.extend(edges_item) + #if faces_item: faces_out.extend(faces_item) self.outputs['Vertices'].sv_set(verts_out)