Dispy benchmark (Delaunay triangulation algorithm)

compute.py

def rgb2gray(rgb):
	r, g, b = rgb[0,:,:], rgb[1,:,:], rgb[2,:,:]
	gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
	return gray


def compute(img):
	try:
#		print('Computing...')
		import time, socket
		host = socket.gethostname()
		import numpy as np
		from scipy import ndimage
		from scipy.ndimage.filters import gaussian_filter
		from scipy.ndimage.measurements import center_of_mass
		from scipy.spatial import Delaunay
		if img.shape != (1024,1024):
			return (host,None,Exception('image is different format',img.shape))
#		from skimage import filter
		min = np.amin(img)
		max = np.amax(img)
		img = gaussian_filter(img,3)
		sz = img.size
		for thres in np.arange(min,max,(max-min)/100):
			im_thresholded = img > thres
			s = np.sum(im_thresholded)
			if s < sz/80:
				break
		im_thresholded = img > thres
		lbl = ndimage.label(im_thresholded)
#		print(lbl)
		num_labels = np.amax(lbl[0])
		ps = center_of_mass(im_thresholded,lbl[0],range(1,num_labels+1))
#		from scipy.misc import toimage
#		toimage(lbl[0]).show()
#		print('%d labels'%(num_labels))
#		print('%d points'%(len(ps)))
		tri = Delaunay(ps)
		import matplotlib.pyplot as plt
#		ps2 = np.array(ps)
#		plt.triplot(ps2[:,0], ps2[:,1], tri.simplices.copy())
#		plt.plot(ps2[:,0], ps2[:,1], 'o', markersize=2)
#		plt.show()
		ds = {}
		from scipy.spatial.distance import cdist
		for i in range(len(tri.simplices)):
			i1 = tri.simplices[i,0]
			i2 = tri.simplices[i,1]
			i3 = tri.simplices[i,2]
			i1,i2,i3 = sorted([i1,i2,i3])
			s1 = '%d-%d'%(i1,i2)
			s2 = '%d-%d'%(i2,i3)
			s3 = '%d-%d'%(i1,i3)
			a = tri.points[i1]
			b = tri.points[i2]
			c = tri.points[i3]
			if(100 < a[0] < 924 and 100 < a[1] < 924):
				if(100 < b[0] < 924 and 100 < b[1] < 924):
					ds[s1] = cdist([a],[b])[0][0]
				if(100 < c[0] < 924 and 100 < c[1] < 924):
					ds[s3] = cdist([a],[c])[0][0]
			if(100 < c[0] < 924 and 100 < c[1] < 924):
				if(100 < b[0] < 924 and 100 < b[1] < 924):
					ds[s2] = cdist([b],[c])[0][0]
		return (host,(tri,ds.values()),None)
	except Exception as e:
		return (host,None,e)


def search():
	import fnmatch
	import os

	matches = []
	for root, dirnames, filenames in os.walk('/Volumes/Macintosh HD/Google Drive/Groves/SEM'):
		for filename in fnmatch.filter(filenames, '*.tiff'):
			matches.append(os.path.join(root, filename))
	return matches


def main(n):
	import time
	import tifffile as tiff
	print('testing compute')
	initiala = time.time()
	files = search()
	files = files[0:n]
	print '%d files.'%(len(files))
	import numpy as np
	spent = []
	for i in range(len(files)):
			try:
				img = tiff.imread(files[i])
				if img.shape != (3,1024,1024):
					continue
				img = rgb2gray(img)
			except:
				continue
			initial = time.time()
			host,res,err = compute(img)
			final = time.time()
	#		print(host,res,err)
			if res is not None:
				spent.append(final - initial)
				ds = res[1]
#				import matplotlib.pyplot as plt
#				plt.hist(ds,60,range=(0,300))
#				plt.show()
#				print('%03d: %d points. %.1f +/- %.1f pixels'%(i,len(res[0].points),np.mean(ds),np.std(ds)))
			else:
				print(err)
	finala = time.time()
	print('Process time: %.2f sec (%d jobs, %.3f sec/job)'%(finala-initiala, len(spent), np.mean(spent)))
	return (finala-initiala, len(spent), np.mean(spent))

if __name__ == '__main__':
	res = []
	for n in [10,30,100,300,1000,3000,10000]:
		res.append(main(n))
	print res