Initial implementation of "Lloyd in Solid" node.

6e806048 · Ilya Portnov · b3d32eb5 · 6e806048 · 6e806048 · 6e806048
--- a/index.md
+++ b/index.md
@@ -321,6 +321,7 @@
    SvVoronoiOnSurfaceNode
    SvLloyd2dNode
    SvLloydOnMeshNode
+    SvLloydSolidNode
    SvRandomPointsOnMesh
    SvFieldRandomProbeMk2Node
    SvPopulateSurfaceNode

--- a/nodes/spatial/lloyd_solid.py
+++ b/nodes/spatial/lloyd_solid.py
+# This file is part of project Sverchok. It's copyrighted by the contributors
+# recorded in the version control history of the file, available from
+# its original location https://github.com/nortikin/sverchok/commit/master
+#  
+# SPDX-License-Identifier: GPL3
+# License-Filename: LICENSE
+
+import bpy
+from bpy.props import FloatProperty, StringProperty, BoolProperty, EnumProperty, IntProperty
+
+from sverchok.node_tree import SverchCustomTreeNode
+from sverchok.data_structure import updateNode, ensure_nesting_level, zip_long_repeat, throttle_and_update_node, get_data_nesting_level
+from sverchok.utils.voronoi3d import lloyd_in_solid
+from sverchok.dependencies import scipy, FreeCAD
+
+if scipy is None or FreeCAD is None:
+    add_dummy('SvLloydSolidNode', "Lloyd in Solid", 'scipy and FreeCAD')
+
+if FreeCAD is not None:
+    import Part
+
+class SvLloydSolidNode(bpy.types.Node, SverchCustomTreeNode):
+    """
+    Triggers: Lloyd Solid
+    Tooltip: Redistribute 3D points in the volume of a Solid body uniformly by use of Lloyd's algorithm
+    """
+    bl_idname = 'SvLloydSolidNode'
+    bl_label = 'Lloyd in Solid'
+    bl_icon = 'OUTLINER_OB_EMPTY'
+    sv_icon = 'SV_VORONOI'
+
+    iterations : IntProperty(
+        name = "Iterations",
+        description = "Number of Lloyd algorithm iterations",
+        min = 0,
+        default = 3,
+        update = updateNode)
+
+    def sv_init(self, context):
+        self.inputs.new('SvSolidSocket', "Solid")
+        self.inputs.new('SvVerticesSocket', "Sites")
+        self.inputs.new('SvStringsSocket', 'Iterations').prop_name = 'iterations'
+        self.outputs.new('SvVerticesSocket', "Sites")
+
+    def process(self):
+
+        if not any(socket.is_linked for socket in self.outputs):
+            return
+
+        solid_in = self.inputs['Solid'].sv_get()
+        sites_in = self.inputs['Sites'].sv_get()
+        iterations_in = self.inputs['Iterations'].sv_get()
+
+        solid_in = ensure_nesting_level(solid_in, 2, data_types=(Part.Shape,))
+        input_level = get_data_nesting_level(sites_in)
+        sites_in = ensure_nesting_level(sites_in, 4)
+        iterations_in = ensure_nesting_level(iterations_in, 2)
+
+        nested_output = input_level > 1
+
+        verts_out = []
+        for params in zip_long_repeat(solid_in, sites_in, iterations_in):
+            new_verts = []
+            for solid, sites, iterations in zip_long_repeat(*params):
+                sites = lloyd_in_solid(solid, sites, iterations)
+                new_verts.append(sites)
+            if nested_output:
+                verts_out.append(new_verts)
+            else:
+                verts_out.extend(new_verts)
+
+        self.outputs['Sites'].sv_set(verts_out)
+
+def register():
+    if scipy is not None and FreeCAD is not None:
+        bpy.utils.register_class(SvLloydSolidNode)
+
+def unregister():
+    if scipy is not None and FreeCAD is not None:
+        bpy.utils.unregister_class(SvLloydSolidNode)
+
--- a/utils/geom.py
+++ b/utils/geom.py
@@ -2165,3 +2165,13 @@ def calc_bounds(vertices, allowance=0):
            y_min - allowance, y_max + allowance,
            z_min - allowance, z_max + allowance)

+TRIVIAL='TRIVIAL'
+def bounding_sphere(vertices, algorithm=TRIVIAL):
+    if algorithm != TRIVIAL:
+        raise Exception("Unsupported algorithm")
+    c = center(vertices)
+    vertices = np.array(vertices) - np.array(c)
+    norms = np.linalg.norm(vertices, axis=1)
+    radius = norms.max()
+    return c, radius
+
--- a/utils/voronoi3d.py
+++ b/utils/voronoi3d.py
@@ -27,11 +27,15 @@ from mathutils.bvhtree import BVHTree
 from sverchok.utils.sv_mesh_utils import mask_vertices, polygons_to_edges
 from sverchok.utils.sv_bmesh_utils import bmesh_from_pydata, pydata_from_bmesh, bmesh_clip
 from sverchok.utils.geom import calc_bounds
-from sverchok.dependencies import scipy
+from sverchok.dependencies import scipy, FreeCAD

 if scipy is not None:
    from scipy.spatial import Voronoi

+if FreeCAD is not None:
+    from FreeCAD import Base
+    import Part
+
 def voronoi3d_layer(n_src_sites, all_sites, make_regions, do_clip, clipping):
    diagram = Voronoi(all_sites)
    src_sites = all_sites[:n_src_sites]
@@ -207,3 +211,48 @@ def lloyd_on_mesh(verts, faces, sites, thickness, n_iterations):

    return points.tolist()

+def lloyd_in_solid(solid, sites, n_iterations, tolerance=1e-4):
+    shell = solid.Shells[0]
+
+    def invert(pt):
+        src = Base.Vector(pt)
+        dist, vs, infos = shell.distToShape(Part.Vertex(src))
+        projection = vs[0][0]
+        dst = src + 2*(projection - src)
+        return (dst.x, dst.y, dst.z)
+        
+    def iteration(pts):
+        n = len(pts)
+        all_pts = pts
+        for pt in pts:
+            if solid.isInside(Base.Vector(pt), tolerance, False):
+                all_pts.append(invert(pt))
+
+        diagram = Voronoi(all_pts)
+        centers = []
+        for site_idx in range(n):
+            region_idx = diagram.point_region[site_idx]
+            region = diagram.regions[region_idx]
+            region_verts = np.array([diagram.vertices[i] for i in region])
+            center = np.mean(region_verts, axis=0)
+            centers.append(tuple(center))
+        return centers
+
+    def restrict(points):
+        result = []
+        for point in points:
+            v = Base.Vector(point)
+            if solid.isInside(v, tolerance, True):
+                result.append(point)
+            else:
+                dist, vs, infos = solid.distToShape(Part.Vertex(v))
+                v = vs[0][0]
+                result.append((v.x, v.y, v.z))
+        return result
+
+    points = restrict(sites)
+    for i in range(n_iterations):
+        points = iteration(points)
+        points = restrict(points)
+    return points
+