Calculating the convex hull of a point data set (Python)

Working with LiDAR point data it was necessary for me to polygonize the point cloud extent. A first approach was to calculate the convex hull of the points. This is predominantly facilitated using scipy spatial’s ConvexHull function.

For my application I required the hull points to be printed out into a txt/csv in order of position (i.e. going clockwise around the hull). This is enabled by simply pulling out the “vertices” from the spatial.ConvexHull object - this gives the indicies of the convex hull points within the original xy input file:

hull = spatial.ConvexHull(xy_file, qhull_options="Qt")
hull_indices = hull.vertices

This is all integrated into a single function:

import os
import sys
import numpy as np
from scipy import spatial

def xy_convex_hull(input_xy_file):
	'''
	Calculates the convex hull of a given xy data set
	returning the indicies of the convex hull points
	in the input data set. A convex hull point
	co-ordinate file is then created using
	write_convex_hull_xy()
	'''
	if os.path.isfile(input_xy_file):
		print "Loading file...."
		xy_file = np.loadtxt(input_xy_file, usecols= (0,1))
	else:
		sys.exit("File for convex hull calculation doesn't exist")

	print "Calculating hull points..."
	hull = spatial.ConvexHull(xy_file, qhull_options="Qt")
	hull_indices = hull.vertices

	print "Found %d hull points" % len(hull_indices)
	print "Hull indicies calculated and now being returned..."

	return xy_file, hull_indices

hull_indicies - which more precisely are the “[i]ndices of points forming the simplical facets of the convex hull” - will look something like this:

array([17218337,    48757,     3476,     1517,      649,       73,
	     38,       37,        1,        0,      209,    27606,
	 153490,   154907,  1422687,  1435840,  1455472,  1933477,
	2183734,  2185942,  2754182,  2977613,  3100991,  7573831,
	8619400, 14809551, 15214362, 16139836, 16146528, 16188976,
       16204639, 16351763, 16649880, 17150945, 17153685, 17260738,
       17268373, 17273051, 17275385, 17299918, 17300751, 17304440,
       17307723, 17318805, 17319327, 17319315, 17319198, 17290692,
       17286814, 17282596, 17279306], dtype=int32)

The hull vertex indices (hull_indices) can then be passed in with the main xy point file (xy_file) (along with an output path (opath) and output file name (file_name) to the below function which writes out the xy positions of hull_indices from xy_file (in the format I required):

def write_convex_hull_xy(xy_file, hull_indices, opath, file_name):
	'''
	Takes the convex hull verticies (an array of indicies) and
	uses them to extract and write out the convex hull vertex
	co-ordinates to a file
	'''

	file_name = "%sconvex_hull_vertices_%s.txt" %(opath, file_name)
	f = open( file_name, 'w')

	for i in range(len(hull_indices)):
		index = hull_indices[i]
		value_y = xy_file[index, 0].astype('float32')
		value_x = xy_file[index, 1].astype('float32')
		f.write("%f %f\n" %(value_y, value_x))

	f.close()

This should give something similar to the following:

540460.812500 7362116.000000
540473.875000 7362116.000000
540473.937500 7362116.000000
540475.687500 7362133.000000
540475.875000 7362134.500000
540478.750000 7362163.000000
	            ... etc.

Ultimately, from this xy point dataset:

Original point image

I managed to get the following:

Convex hull point verticies

Written on November 6, 2014