Constructing
3D Models
from
Voronoi Diagrams

 




Introduction

Qhull is a program to compute convex hulls, Delaunay triangulations, Voronoi diagrams, halfspace intersections about a point, furthest-site Delaunay triangulations, and furthest-site Voronoi diagrams. Its output can be used to create 3D models from the Voronoi cells.

The program rbox, which is part of the Qhull package, can generate random or regular points (also called generators). When rbox is invoked as follows:

C:\>rbox 10 D2

it will generate 10 random 2D points in the cube (-0.5,-0.5)-(0.5,0.5):

2 rbox 10 D2
10
-0.02222276248244826 -0.4979727817680433
-0.4285431913366012 0.4745826469497594
0.3105396575392593 0.2400179190933871
-0.01883958887200765 0.3630260628303755
0.3790312361708201 0.3779794437605696
-0.2994955874043476 0.3776609263174803
0.3471817493878135 0.08365533089605659
-0.00485819764887746 0.3482682405489201
0.3443122672329771 -0.1437312230875075
0.309330780347186 -0.07758103877080702

When rbox is invoked as follows:

C:\>rbox 10

it will generate 10 random 3D points in the unit cube centered at the origin:

3 rbox 10
10
-0.374949442106252 0.2247255418679914 -0.03782349642163468
0.3004920247946792 0.3694544596313075 0.4210962191420572
0.3641479121187217 0.2339522165562512 0.03489791698278666
-0.4707134565997063 -0.2810653003687648 0.1354949890035158
0.2642752088273644 -0.3265712199048803 0.3175057017407433
0.3183227589524563 0.05060330969337679 0.4898217078203537
0.4334355904100794 -0.2480392830893763 0.2037671452423252
-0.2855954200826543 -0.002227006947833121 -0.4293096679535784
-0.4075898489985521 -0.3625928418362446 -0.09789381697577781
-0.3013850588364387 -0.3786854184034145 0.4341719447953365

Constructing Meshes from 2D Voronoi Diagrams

The following image shows a 2D Voronoi diagram as a 3D mesh:

To create the mesh, we have to do to the following. First, we generate 1000 random 2D points. The output of rbox is sent immediately to qvoronoi. In turn, the output of qvoronoi is written to the file out.txt:

C:\>rbox 1000 D2|qvoronoi p FN >out.txt

The option ‘p’ prints a list of Voronoi vertices. The option ‘FN’ outputs the count and Voronoi vertices — as indices — for each Voronoi region. The contents of the file out.txt looks as follows:

2
1986
-0.1567880435769815 -75.03601498115664
-0.06742416186229336 -44.60355654095014
0.4576158688630084 -1.152698281169307
0.03696212219638438 1.341102309313179
-0.02692708739876709 3.667325402516894
0.09446768016164753 21.99531430484879
0.4961194581483074 0.02083787067511567
-0.6501734950374665 -0.4344678811656289
.
.
.
1000
11 1383 83 271 -2130 7 888 889 453 79 80 1382
6 1284 182 336 335 789 1282
5 360 145 56 55 359
4 1040 1038 1037 1039
5 1795 1098 1102 1101 1793
4 1961 1959 1958 1960
5 196 67 65 60 195
6 546 249 252 251 250 545
6 675 138 41 12 10 674
.
.
.

There are 1986 vertices and 1000 cells. Each region line contains the count of vertices (e.g. 11 for the first region) and the vertex indices. Regions with negative indices will be ignored.

We can now use the following Python script to convert the Voronoi diagram to a 3D mesh. The output of the script is a Wavefront obj file.

def CrossProduct(a, b):
    x = (a[1] * b[2]) - (a[2] * b[1])
    y = -((a[2] * b[0]) - (a[0] * b[2]))
    z = (a[0] * b[1]) - (a[1] * b[0])
    return [x, y, z]

def IsBounded(item):
    for x in item:
        if x < 0:
            return False
    return True
    
class VoronoiDiagram:
    def __init__(self):
        self.Vertices = []
        self.Regions = []
    def UpdateVertexOrder(self):
        for r in filter(IsBounded, self.Regions):
            v0 = self.Vertices[r[0]]
            v1 = self.Vertices[r[1]]
            v2 = self.Vertices[r[2]]
            a = (v0[0] - v1[0], 0, v0[1] - v1[1])
            b = (v2[0] - v1[0], 0, v2[1] - v1[1])
            t = CrossProduct(a, b)
            if t[1] < 0:
                r.reverse()
    def Load(self, filename):
        del self.Vertices[:]
        del self.Regions[:]
        f = open(filename, 'r')
        f.readline()
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            x = float(t[0])
            y = float(t[1])
            self.Vertices.append([x, y])
            i = i - 1
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            n = int(t[0])
            region = []
            j = 0
            while j < n:
                region.append(int(t[1 + j]))
                j = j + 1
            self.Regions.append(region)
            i = i - 1
        f.close()
        self.UpdateVertexOrder()
    def SaveAsObj(self, filename):
        f = open(filename, 'w')
        for v in self.Vertices:
            s = 'v ' + str(v[0]) + ' 0 ' + str(v[1]) + '\n'
            f.write(s)
        for region in filter(IsBounded, self.Regions):
            s = 'f '
            for v in region:
                s += str(v + 1) + ' '
            s += '\n'
            f.write(s)
        f.close()

v = VoronoiDiagram()
v.Load('out.txt')
v.SaveAsObj('obj.obj')

This script does little more than parsing the file out.txt and saving the Voronoi regions as polygons. All vertices will have a y-coordinate of 0. All normals will point straight upwards because we change the vertex order if necessary (because we cannot trust the vertex order outputted by qhull, we compute the cross product of two vectors of a polygon; if the y-coordinate is negative, we reverse the order of the vertex indices). If you load the file obj.obj in your favorite 3D graphics program, you should see something like the image above. IF YOU DON'T SEE ANYTHING, OR THE MODEL APPEARS VERY SMALL, SCALE YOUR MODEL! This is perfectly normal because all coordinates are in the range [-0.5, 0.5].

Why are the outer regions so spiky? Let's take a look at a Voronoi diagram:

As you can see, some of the outer regions have edges that converge, others have edges that diverge. The regions with diverging edges are unbounded and aren't even included in the mesh. The regions with converging edges are included, but they are not clipped by a bounding box!. A good 3D graphics program should be able to clip the mesh.

Constructing Polyhedrons from 3D Voronoi Diagrams

Given a set of 3D points, how can we construct Voronoi polyhedrons so that we can view them as 3D solid models?

First, we're going to create a file with 1000 random points:

C:\>rbox 1000 >points.txt

The output of rbox is redirected to the file points.txt and looks as follows:

3 rbox 1000
1000
-0.3644684295770418 0.3791030383474222 -0.41524137502093
0.03820936012846354 0.1847114671810637 0.4456235537730376
-0.4049391373106643 0.1879184760040775 0.3458208170214862
0.2104650607430051 0.2862703476904616 0.3457274803386325
-0.358244567511831 -0.01644728054892952 -0.4294479702873602
0.2679628285281014 -0.3487469384900731 -0.3897963863702457
-0.3078665871246407 -0.3137313074560196 0.1179141291584019
-0.2172360706303604 -0.08664129729069892 -0.1802867997254122
-0.08024548700102185 0.3140966890511072 0.02304651078120479
0.3427026065464156 -0.1972983695541493 0.00630053459322133
.
.
.

Now we're going to calculate the vertices and faces:

C:\>type points.txt|qvoronoi p Fv >out.txt

The contents of the file out.txt looks as follows:

3
6298
7.851314493941277 -0.1480177141824365 -0.6230659162387991
0.442558342462492 0.4568540929808236 0.1480251963386509
0.07575807322885306 0.4463396644341515 0.09236281428869744
-1.265510019293124 0.2299237088889936 -0.4370463335702945
-0.6613388325413993 0.3848427453253555 -0.344975621925081
-1.137572220612269 0.2830507041353281 -0.2742760853332658
-1.532407384758669 0.3227718805850115 -0.3564498198927059
-2.404008618045335 0.3646814753565159 -0.252538955426773
0.235578274785728 0.9923899276802044 0.1729807422096146
0.3581420751492493 5.04987634278121 0.558376131829994
.
.
.
7352
9 0 393 756 1601 2848 2849 763 758 762
5 0 584 756 760 762
10 0 712 756 760 761 757 5160 5161 5126 1601
6 0 942 757 5160 5158 1208
7 0 703 757 761 759 1207 1208
8 0 738 758 759 1207 2845 2844 763
7 0 855 758 762 760 761 759
5 0 784 763 2849 2844
9 0 67 1207 2845 2851 2885 5162 5158 1208
8 0 575 1601 5126 5127 2885 2851 2848
.
.
.

The number 6298 on line 2 is the number of vertices. It is obvious that each vertex has a x, y and z-coordinate. The number 7352 is the number of faces. Let's take a look at the first face. The number 9 is the count of indices. The second pair of indices, 0 and 393, indcates a generator pair: 0 is the first point in the file points.txt and 393 is the 394th point (if we start counting from 1 instead of 0). The next seven indices are the indices of the vertices of the face: e.g. 756 is the 756th vertex (not in the file points.txt of course, but in the list of vertices in the file out.txt). An index of 0 indicates the vertex-at-infinity (we're going to ignore polyhedrons which have faces with a vertex-at-infinity).

Armed with these two files (points.txt and out.txt), we are ready to create solid 3D models. The following Python script will create a Wavefront obj file from our Voronoi data:

def GetPlaneEquation(v1, v2, v3):
    x1, y1, z1 = v1[0], v1[1], v1[2]
    x2, y2, z2 = v2[0], v2[1], v2[2]
    x3, y3, z3 = v3[0], v3[1], v3[2]
    A = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
    B = z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2)
    C = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    D = -(x1 * (y2 * z3 - y3 * z2) + \
        x2 * (y3 * z1 - y1 * z3) + \
        x3 * (y1 * z2 - y2 * z1))
    return [A, B, C, D]
    
class VoronoiDiagram:
    def __init__(self):
        self.Generators = []
        self.Vertices = []
        self.Faces = []
    def Load(self, pointsfilename, qvfilename):
        del self.Generators[:]
        f = open(pointsfilename, 'r')
        l = f.readline()
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            x = float(t[0])
            y = float(t[1])
            z = float(t[2])
            self.Generators.append([x, y, z])
            i = i - 1
        f.close()
        del self.Vertices[:]
        del self.Faces[:]
        f = open(qvfilename, 'r')
        l = f.readline()
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            x = float(t[0])
            y = float(t[1])
            z = float(t[2])
            self.Vertices.append([x, y, z])
            i = i - 1
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            n = int(t[0])
            face = []
            j = 0
            while j < n:
                face.append(int(t[1 + j]))
                j = j + 1
            self.Faces.append(face)
            i = i - 1
        f.close()
    def GetFacesForGenerator(self, x):
        faces = []
        i = 0
        while i < len(self.Faces):
            if self.Faces[i][0] == x or self.Faces[i][1] == x:
                faces.append(self.Faces[i][2:])
            i = i + 1
        return faces
    def HasVertexAtInfinity(self, faces):
        for face in faces:
            for index in face:
                if index == 0:
                    return True
        return False
    def UpdateNormal(self, g, face):
        v1 = self.Vertices[face[0] - 1]
        v2 = self.Vertices[face[1] - 1]
        v3 = self.Vertices[face[2] - 1]
        A, B, C, D = GetPlaneEquation(v1, v2, v3)
        x, y, z = g
        if A * x + B * y + C * z + D > 0:
            face.reverse()
    def SaveAsObj(self, filename):
        f = open(filename, 'w')
        base = 1
        i = 0
        while i < len(self.Generators):
            faces = self.GetFacesForGenerator(i)
            if not self.HasVertexAtInfinity(faces):
                f.write('g cell%d\n' % i)
                vertices = []
                for face in faces:
                    for index in face:
                        if vertices.count(index) == 0:
                            vertices.append(index)
                for index in vertices:
                    f.write('v %.12f %.12f %.12f\n' % \
                            tuple(self.Vertices[index - 1]))
                for face in faces:
                    self.UpdateNormal(self.Generators[i], face)
                    s = 'f '
                    for index in face:
                        s += str(base + (vertices.index(index))) + ' '
                    s += '\n'
                    f.write(s)
                f.write('\n')
                base = base + len(vertices)
            i = i + 1
        f.close()

v = VoronoiDiagram()
v.Load('points.txt', 'out.txt')
v.SaveAsObj('obj.obj')

This Python script defines a class VoronoiDiagram which can parse the files points.txt and out.txt and save the Voronoi diagram as an .obj file. The function GetFacesForGenerator() looks at the list of faces and returns the faces (only the vertex indices) where the generator pair includes the given generator. The function HasVertexAtInfinity() will return True if there's a face in the list of faces which has a vertex index of 0 (the vertex-at-infinity). Polyhedrons with such faces will not be converted to a 3D solid model. The function SaveAsObj() iterates over the generators. If the polyhedron of the generator doesn't have a vertex-at-infinity, the function creates a list of vertices (a separate list per polyhedron is created so that the polyhedrons do not share their vertices and faces; this way, they can be move and sized independently). Once the list of vertices is constructed, the faces can be written. The normal of the face is updated if necessary: if the generator lies at the same side as the normal, the normal is flipped (in the Voronoi data in the file out.txt, the faces are shared between generators, so for one generator, the normal will be ok, but for the adjacent generator, the normal will be incorrect because the indices of the vertices will be in clockwise order).

Here is a wire frame model of the Voronoi diagram:

And here is a solid model:

Again, the outer regions are spiky. See Constructing Meshes from 2D Voronoi Diagrams for an explanation.

Constructing Manifolds from 3D Voronoi Diagrams

If you want to remove the invisible, inner polygons so that there is only one polyhedron, use the following Python script:

def GetPlaneEquation(v1, v2, v3):
    x1, y1, z1 = v1[0], v1[1], v1[2]
    x2, y2, z2 = v2[0], v2[1], v2[2]
    x3, y3, z3 = v3[0], v3[1], v3[2]
    A = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
    B = z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2)
    C = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    D = -(x1 * (y2 * z3 - y3 * z2) + \
        x2 * (y3 * z1 - y1 * z3) + \
        x3 * (y1 * z2 - y2 * z1))
    return [A, B, C, D]
    
class VoronoiDiagram:
    def __init__(self):
        self.Generators = []
        self.Vertices = []
        self.Faces = []
    def Load(self, pointsfilename, qvfilename):
        del self.Generators[:]
        f = open(pointsfilename, 'r')
        l = f.readline()
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            x = float(t[0])
            y = float(t[1])
            z = float(t[2])
            self.Generators.append([x, y, z])
            i = i - 1
        f.close()
        del self.Vertices[:]
        del self.Faces[:]
        f = open(qvfilename, 'r')
        l = f.readline()
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            x = float(t[0])
            y = float(t[1])
            z = float(t[2])
            self.Vertices.append([x, y, z])
            i = i - 1
        l = f.readline()
        i = int(l.strip())
        while i > 0:
            l = f.readline()
            t = l.split()
            n = int(t[0])
            face = []
            j = 0
            while j < n:
                face.append(int(t[1 + j]))
                j = j + 1
            self.Faces.append(face)
            i = i - 1
        f.close()
    def GetFacesForGenerator(self, x):
        faces = []
        i = 0
        while i < len(self.Faces):
            if self.Faces[i][0] == x or self.Faces[i][1] == x:
                faces.append(self.Faces[i][2:])
            i = i + 1
        return faces
    def HasVertexAtInfinity(self, faces):
        for face in faces:
            for index in face:
                if index == 0:
                    return True
        return False
    def UpdateNormal(self, g, face):
        v1 = self.Vertices[face[0] - 1]
        v2 = self.Vertices[face[1] - 1]
        v3 = self.Vertices[face[2] - 1]
        A, B, C, D = GetPlaneEquation(v1, v2, v3)
        x, y, z = g
        if A * x + B * y + C * z + D > 0:
            face.reverse()
    def IsFaceVisible(self, face):
        faces0 = self.GetFacesForGenerator(face[0])
        faces1 = self.GetFacesForGenerator(face[1])
        b0 = self.HasVertexAtInfinity(faces0)
        b1 = self.HasVertexAtInfinity(faces1)
        if b0 and not b1:
            return 1
        elif b1 and not b0:
            return 0
        return -1
    def SaveVertexIndices(self, vilist, indices):
        for index in indices:
            if vilist.count(index) == 0:
                vilist.append(index)
    def WriteObjVertices(self, f, vilist):
        for index in vilist:
            f.write('v %.12f %.12f %.12f\n' % tuple(self.Vertices[index - 1]))
    def WriteObjFaces(self, f, vilist, faces):
        for face in faces:
            f.write('f ')
            for index in face:
                f.write(str(vilist.index(index) + 1) + ' ')
            f.write('\n')
    def SaveAsObj(self, filename):
        objvertexindices = []
        objfaces = []
        for face in self.Faces:
            n = self.IsFaceVisible(face)
            if n != -1:
                indices = face[2:]
                if not self.HasVertexAtInfinity((indices,)):
                    self.UpdateNormal(self.Generators[face[n]], indices)
                    self.SaveVertexIndices(objvertexindices, indices)
                    objfaces.append(indices)
        f = open(filename, 'w')
        self.WriteObjVertices(f, objvertexindices)
        self.WriteObjFaces(f, objvertexindices, objfaces)
        f.close()

v = VoronoiDiagram()
v.Load('points.txt', 'out.txt')
v.SaveAsObj('obj.obj')

Here is an image of the 3D model:

Downloads

The file voronoi3d.zip contains the Python scripts convert2d.py, convert3d.py and manifold3d.py, and three 3D models: obj2d.obj, obj3d.obj and objm3d.obj.

 





 

Back