Tuesday, March 1, 2011

Rough Draft: How to generate 2D Delaunay Triangulations using Qhull and Python

Delaunay triangulations partition a space into regions. Delaunay triangulations can be useful for interpolation and visualizing spatial data.

sample_delaunay

Qhull is a program which can generate tesselations, convex hulls, and vonoroi diagrams from a set of points. This program is available as a precompiled executable and source code. By interfacing to the command line version of this program, a Delaunay triangulation can be generated.

Sample CodE

To use this code, download Qhull and copy the ‘qhull.exe’ to the working directory. Make sure the code has read, write, and execute privileges in the working directory. Make sure there are not files named ‘data.txt’ or ‘results.txt’ which need to be preserved.

'''
This module generates 2D tesselations for a set of points.
The vornoi cells can be filtered by supplying a value associated with 
   each node.
'''

__author__ = 'Ed Tate'
__email__  = 'edtate-at-gmail-dot-com'

def delaunay2D(xpt,ypt,cpt=None,threshold=0):
    
    if cpt is None:
        cpt = [-1 for x in xpt]    
        
    # write the data file
    pts_filename = 'data.txt'
    pts_F = open(pts_filename,'w')
    #print pts_F
    pts_F.write('2 # this is a 2-D input set\n')
    pts_F.write('%i # number of points\n' % len(xpt))
    for i,(x,y) in enumerate(zip(xpt,ypt)):
        pts_F.write('%f %f # data point %i\n' % (x,y,i))
    pts_F.close()

    # trigger the shell command
    import subprocess

    p = subprocess.Popen('qhull TI data.txt TO results.txt d i Qc Qt Qbb', shell=True)
    p.wait()

    # open the results file and parse results
    results = open('results.txt','r')
    print results

    # get 'i' results
    data = results.readline()
    tri_list = []
    numLines = int(data)
    print numLines
    for i in range(numLines):
        # load each triplet of indexes to x,y points
        data = results.readline()
        idx1,idx2,idx3,dummy = data.split(' ')
        idx1 = int(idx1)
        idx2 = int(idx2)
        idx3 = int(idx3)
        tri_list.append([idx1,idx2,idx3])
        
    #################
    #this generates a fillable collection of tesselations
    x_list = []
    y_list = []   
    
    for t in tri_list:
        if all([ cpt[t[i]]<threshold for i in range(3)]):
            short_x_list = [ xpt[t[i]] for i in range(3)]
            short_y_list = [ ypt[t[i]] for i in range(3)]
            x_list.extend(short_x_list)
            x_list.append(xpt[t[0]])
            x_list.append(None)
            y_list.extend(short_y_list)
            y_list.append(ypt[t[0]])
            y_list.append(None)
    
    return (x_list,y_list)


if __name__=='__main__':
    
    import random
    
    N=100
    xpt = [random.random()-0.5 for i in range(0,N)]
    ypt = [random.random()-0.5 for i in range(0,N)]
    cpt = [random.random()-0.5 for i in range(0,N)]
    
    import matplotlib.pyplot as pp
    
    pp.figure()
    (x_list,y_list) = delaunay2D(xpt,ypt,cpt)
    pp.plot(xpt,ypt,'b.')
    pp.plot(x_list,y_list,'k-')

    pp.fill(x_list,y_list,'g',alpha=0.25,edgecolor='none')    

    pp.show()
    

When run, this module will produce a diagram with the cells with filled triangulation where the vertices have a random value greater than 0.

sample_delaunay_2

Test Conditions


Answers

  • How to tesselate a set of points
  • How to generate a Delaunay Triangulation
  • How to use Qhull with Python
This work is licensed under a Creative Commons Attribution By license.

No comments:

Post a Comment