Skip to content
Snippets Groups Projects
delaunay_voronoi.py 11.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • # -*- coding:utf-8 -*-
    
    
    # ##### BEGIN GPL LICENSE BLOCK #####
    #
    #  This program is free software; you can redistribute it and/or
    #  modify it under the terms of the GNU General Public License
    #  as published by the Free Software Foundation; either version 2
    #  of the License, or (at your option) any later version.
    #
    #  This program is distributed in the hope that it will be useful,
    #  but WITHOUT ANY WARRANTY; without even the implied warranty of
    #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    #  GNU General Public License for more details.
    #
    #  You should have received a copy of the GNU General Public License
    #  along with this program; if not, write to the Free Software Foundation,
    #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
    #
    # ##### END GPL LICENSE BLOCK #####
    
    bl_info = {
        "name": "Delaunay Voronoi",
        "description": "Points cloud Delaunay triangulation in 2.5D "
                       "(suitable for terrain modelling) or Voronoi diagram in 2D",
        "author": "Domlysz, Oscurart",
        "version": (1, 3),
        "blender": (2, 7, 0),
    
        "location": "3D View > Toolshelf > Create > Delaunay Voronoi",
    
        "warning": "",
        "wiki_url": "https://github.com/domlysz/BlenderGIS/wiki",
        "category": "Add Mesh"
        }
    
    
    import bpy
    from .DelaunayVoronoi import (
    
            computeVoronoiDiagram,
            computeDelaunayTriangulation,
            )
    
    from bpy.types import (
            Operator,
            Panel,
            )
    from bpy.props import EnumProperty
    
    
    try:
        from scipy.spatial import Delaunay
        import bmesh
        import numpy as np
        HAS_SCIPY = True
    except:
        HAS_SCIPY = False
        pass
    
    
    # Globals
    # set to True to enable debug_prints
    DEBUG = False
    
    
    def debug_prints(text=""):
        global DEBUG
        if DEBUG and text:
            print(text)
    
    
    
    class Point:
        def __init__(self, x, y, z):
            self.x, self.y, self.z = x, y, z
    
    
    def unique(L):
        """Return a list of unhashable elements in s, but without duplicates.
        [[1, 2], [2, 3], [1, 2]] >>> [[1, 2], [2, 3]]"""
        # For unhashable objects, you can sort the sequence and
        # then scan from the end of the list, deleting duplicates as you go
        nDupli = 0
        nZcolinear = 0
        # sort() brings the equal elements together; then duplicates
        # are easy to weed out in a single pass
        L.sort()
        last = L[-1]
        for i in range(len(L) - 2, -1, -1):
            if last[:2] == L[i][:2]:    # XY coordinates compararison
                if last[2] == L[i][2]:  # Z coordinates compararison
                    nDupli += 1         # duplicates vertices
                else:  # Z colinear
                    nZcolinear += 1
                del L[i]
            else:
                last = L[i]
        # list data type is mutable, input list will automatically update
        # and doesn't need to be returned
        return (nDupli, nZcolinear)
    
    
    def checkEqual(lst):
        return lst[1:] == lst[:-1]
    
    
    class ToolsPanelDelaunay(Panel):
        bl_category = "Create"
        bl_label = "Delaunay Voronoi"
        bl_space_type = "VIEW_3D"
        bl_context = "objectmode"
        bl_region_type = "TOOLS"
        bl_options = {"DEFAULT_CLOSED"}
    
        def draw(self, context):
            layout = self.layout
    
    
            box = layout.box()
            col = box.column(align=True)
            col.label("Point Cloud:")
            col.operator("delaunay.triangulation")
            col.operator("voronoi.tesselation")
    
    
    
    class OBJECT_OT_TriangulateButton(Operator):
        bl_idname = "delaunay.triangulation"
        bl_label = "Triangulation"
    
        bl_description = ("Terrain points cloud Delaunay triangulation in 2.5D\n"
                          "Needs an existing Active Mesh Object")
        bl_options = {"REGISTER", "UNDO"}
    
        @classmethod
        def poll(cls, context):
            obj = context.active_object
            return (obj is not None and obj.type == "MESH")
    
        def execute(self, context):
            # move the check into the poll
            obj = context.active_object
    
            if HAS_SCIPY:
                # Use scipy when present (~18 x faster)
                bpy.ops.object.mode_set(mode='EDIT')
                bm = bmesh.from_edit_mesh(obj.data)
                points_3D = [list(v.co) for v in bm.verts]
                points_2D = np.array([[v[0], v[1]] for v in points_3D])
                print("Triangulate " + str(len(points_3D)) + " points...")
                # Triangulate
                tri = Delaunay(points_2D)
                faces = tri.simplices.tolist()
                # Create new mesh structure
                print("Create mesh...")
                bpy.ops.object.mode_set(mode='OBJECT')
                mesh = bpy.data.meshes.new("TIN")
                mesh.from_pydata(points_3D, [], faces)
                mesh.update(calc_edges=True)
                my = bpy.data.objects.new("TIN", mesh)
                context.scene.objects.link(my)
                my.matrix_world = obj.matrix_world.copy()
                obj.select = False
                my.select = True
                context.scene.objects.active = my
                self.report({'INFO'}, "Mesh created (" + str(len(faces)) + " triangles)")
                print("Total :%s faces  %s verts" % (len(faces), len(points_3D)))
                return {'FINISHED'}
    
    
            # Get points coodinates
            r = obj.rotation_euler
            s = obj.scale
            mesh = obj.data
            vertsPts = [vertex.co for vertex in mesh.vertices]
    
            # Remove duplicate
            verts = [[vert.x, vert.y, vert.z] for vert in vertsPts]
            nDupli, nZcolinear = unique(verts)
            nVerts = len(verts)
    
    
            debug_prints(text=str(nDupli) + " duplicate points ignored")
            debug_prints(str(nZcolinear) + " z colinear points excluded")
    
    
                self.report({"WARNING"},
                            "Not enough points to continue. Operation Cancelled")
    
                return {"CANCELLED"}
    
    
            # Check colinear
            xValues = [pt[0] for pt in verts]
            yValues = [pt[1] for pt in verts]
    
            if checkEqual(xValues) or checkEqual(yValues):
                self.report({'ERROR'}, "Points are colinear")
                return {'FINISHED'}
    
            # Triangulate
    
            debug_prints(text="Triangulate " + str(nVerts) + " points...")
    
    
            vertsPts = [Point(vert[0], vert[1], vert[2]) for vert in verts]
            triangles = computeDelaunayTriangulation(vertsPts)
            # reverse point order --> if all triangles are specified anticlockwise then all faces up
            triangles = [tuple(reversed(tri)) for tri in triangles]
    
    
            debug_prints(text=str(len(triangles)) + " triangles")
    
    
            # Create new mesh structure
    
            debug_prints(text="Create mesh...")
    
            tinMesh = bpy.data.meshes.new("TIN")        # create a new mesh
            tinMesh.from_pydata(verts, [], triangles)   # Fill the mesh with triangles
            tinMesh.update(calc_edges=True)             # Update mesh with new data
    
            # Create an object with that mesh
            tinObj = bpy.data.objects.new("TIN", tinMesh)
    
            # Place object
            tinObj.location = obj.location.copy()
            tinObj.rotation_euler = r
            tinObj.scale = s
    
            # Update scene
            bpy.context.scene.objects.link(tinObj)  # Link object to scene
            bpy.context.scene.objects.active = tinObj
            tinObj.select = True
            obj.select = False
    
    
            self.report({"INFO"},
                         "Mesh created (" + str(len(triangles)) + " triangles)")
    
    
            return {'FINISHED'}
    
    
    class OBJECT_OT_VoronoiButton(Operator):
        bl_idname = "voronoi.tesselation"
        bl_label = "Diagram"
    
        bl_description = ("Points cloud Voronoi diagram in 2D\n"
                          "Needs an existing Active Mesh Object")
    
        bl_options = {"REGISTER", "UNDO"}
    
        meshType = EnumProperty(
    
                items=[('Edges', "Edges", "Edges Only - do not fill Faces"),
                       ('Faces', "Faces", "Fill Faces in the new Object")],
                name="Mesh type",
                description="Type of geometry to generate"
                )
    
        @classmethod
        def poll(cls, context):
            obj = context.active_object
            return (obj is not None and obj.type == "MESH")
    
        def execute(self, context):
            # move the check into the poll
            obj = context.active_object
    
    
            # Get points coodinates
            r = obj.rotation_euler
            s = obj.scale
            mesh = obj.data
            vertsPts = [vertex.co for vertex in mesh.vertices]
    
            # Remove duplicate
            verts = [[vert.x, vert.y, vert.z] for vert in vertsPts]
            nDupli, nZcolinear = unique(verts)
            nVerts = len(verts)
    
    
            debug_prints(text=str(nDupli) + " duplicates points ignored")
            debug_prints(text=str(nZcolinear) + " z colinear points excluded")
    
                self.report({"WARNING"},
                            "Not enough points to continue. Operation Cancelled")
    
                return {"CANCELLED"}
    
    
            # Check colinear
            xValues = [pt[0] for pt in verts]
            yValues = [pt[1] for pt in verts]
    
            if checkEqual(xValues) or checkEqual(yValues):
    
                self.report({"WARNING"},
                            "Points are colinear. Operation Cancelled")
    
                return {"CANCELLED"}
    
            debug_prints(text="Tesselation... (" + str(nVerts) + " points)")
    
    
            xbuff, ybuff = 5, 5
            zPosition = 0
            vertsPts = [Point(vert[0], vert[1], vert[2]) for vert in verts]
    
            if self.meshType == "Edges":
                pts, edgesIdx = computeVoronoiDiagram(
    
                                    vertsPts, xbuff, ybuff,
                                    polygonsOutput=False, formatOutput=True
    
                                    )
            else:
                pts, polyIdx = computeVoronoiDiagram(
                                    vertsPts, xbuff, ybuff, polygonsOutput=True,
                                    formatOutput=True, closePoly=False
                                    )
    
            pts = [[pt[0], pt[1], zPosition] for pt in pts]
    
            # Create new mesh structure
            voronoiDiagram = bpy.data.meshes.new("VoronoiDiagram")  # create a new mesh
    
            if self.meshType == "Edges":
                # Fill the mesh with triangles
                voronoiDiagram.from_pydata(pts, edgesIdx, [])
            else:
                # Fill the mesh with triangles
                voronoiDiagram.from_pydata(pts, [], list(polyIdx.values()))
    
            voronoiDiagram.update(calc_edges=True)  # Update mesh with new data
            # create an object with that mesh
            voronoiObj = bpy.data.objects.new("VoronoiDiagram", voronoiDiagram)
            # place object
            voronoiObj.location = obj.location.copy()
            voronoiObj.rotation_euler = r
            voronoiObj.scale = s
    
            # update scene
            bpy.context.scene.objects.link(voronoiObj)  # Link object to scene
            bpy.context.scene.objects.active = voronoiObj
            voronoiObj.select = True
            obj.select = False
    
            # Report
            if self.meshType == "Edges":
    
                self.report({"INFO"}, "Mesh created (" + str(len(edgesIdx)) + " edges)")
    
                self.report({"INFO"}, "Mesh created (" + str(len(polyIdx)) + " polygons)")
    
    # Register
    def register():
        bpy.utils.register_class(OBJECT_OT_VoronoiButton)
        bpy.utils.register_class(OBJECT_OT_TriangulateButton)
        bpy.utils.register_class(ToolsPanelDelaunay)
    
    
    def unregister():
        bpy.utils.unregister_class(OBJECT_OT_VoronoiButton)
        bpy.utils.unregister_class(OBJECT_OT_TriangulateButton)
        bpy.utils.unregister_class(ToolsPanelDelaunay)
    
    
    if __name__ == "__main__":
        register()