--- a
+++ b/10x/batch.py
@@ -0,0 +1,309 @@
+import os
+import sys
+import pickle
+
+import numpy as np
+import matplotlib
+
+import skimage as ski
+from skimage import io, filter, color, exposure, morphology, feature, draw, measure, transform
+
+from scipy.spatial import distance
+from scipy import ndimage
+
+
+
+
+
+
+
+# Things to measure for each image
+in_count = []
+out_count = []
+in_size = []
+out_size = []
+in_area_density = []
+out_area_density = []
+in_point_density = []
+out_point_density = []
+var_point_density = []
+cov_point_density = []
+var_area_density = []
+cov_area_density = []
+fs_area = []
+pairwise = []
+name = []
+imageid = []
+
+
+
+
+
+
+
+F = open("filelist%i.p" % int(sys.argv[1]), "r")
+files = pickle.load(F)
+F.close()
+
+
+
+
+
+
+
+
+for f in files :
+
+	print "Thread %s, image %s" % (sys.argv[1], f)
+
+    ## ID
+	name.append(int(f.split("data/Sheep")[1].split("-")[0].split(".jpg")[0]))
+	imageid.append(int(f.split("data/Sheep")[1].split("-")[2].split(".jpg")[0]))
+
+
+
+
+
+
+
+	# Process for nuclei
+
+	A = io.imread(f)
+	print "Loaded image."
+	
+	binary = filter.threshold_adaptive(exposure.adjust_sigmoid(A[:, :, 0], cutoff=0.4, gain = 30), 301).astype(bool)
+	print "Binarised."
+
+	clean = morphology.binary_closing(binary, morphology.disk(3)).astype(bool)
+	clean = morphology.remove_small_objects(clean, 200)
+	clean = morphology.remove_small_objects( (1-clean).astype(bool), 200)
+	print "Cleaned."
+
+
+
+
+
+
+
+
+
+	# Find contour of inflammatory zone
+
+	local_density = filter.gaussian_filter(clean, 61)
+	local_density -= local_density.min()
+	local_density /= local_density.max()
+
+	ent = filter.gaussian_filter(filter.rank.entropy(local_density, morphology.disk(3)), 75)
+	ent -= ent.min()
+	ent /= ent.max()
+
+	info = ent * (1 + local_density)
+
+	bw = (info) > filter.threshold_otsu(info)
+
+	C = measure.find_contours(bw, 0.5)
+	centroid = []
+	vals = []
+
+	for c in C :
+	    centroid.append(np.linalg.norm([c[:, 1].mean() - bw.shape[1] / 2, c[:, 0].mean() - bw.shape[0] / 2]))
+	    vals.append(local_density.T[c.astype(int)].sum())
+
+	cent = C[np.argmin(centroid / np.array(vals))]
+	contour = matplotlib.path.Path(cent)
+
+	print "Contour of inflammatory zone found."
+
+
+
+
+
+
+
+
+	# Segmentation
+
+	dist = ndimage.distance_transform_edt(clean)
+	peaks = feature.peak_local_max(dist, indices=False, labels = clean)
+	markers = ndimage.label(peaks)[0]
+	labels = morphology.watershed(-dist, markers, mask=clean)
+	print "Image segmented."
+
+
+
+
+
+
+
+
+
+
+
+	# Properties of each nucleus
+	nuclei = measure.regionprops(labels)
+	contour = matplotlib.path.Path(cent, closed=True)
+
+	incount = 0
+	outcount = 0
+	insize = []
+	outsize = []
+	outer_nuclei = []
+
+	# Number of nuclei, and their sizes, inside and outside the focus
+	for n in nuclei :
+	    if contour.contains_point(n.centroid) :
+	        incount += 1
+	        insize.append((labels == n.label).sum())
+	    else :
+	        outcount += 1
+	        outsize.append((labels == n.label).sum())
+	        outer_nuclei.append(n)
+
+	in_count.append(incount)
+	out_count.append(outcount)
+	in_size.append(np.mean(insize))
+	out_size.append(np.mean(outsize))
+	print "Nuclei identified."
+
+
+
+
+
+
+
+
+
+
+
+
+	# Using the shoelace algorithm, compute the area of the focus
+	x = contour.vertices[:, 0]
+	y = contour.vertices[:, 1]
+	x = np.append(x, contour.vertices[0, 0])
+	y = np.append(y, contour.vertices[0, 1])
+	fsarea = np.abs(np.sum(x[:-1] * y[1:] - x[1:] * y[:-1])) / (2. * 3456 * 5184)
+	fs_area.append(fsarea)
+
+
+
+
+
+
+
+
+
+
+
+
+
+	# Area densities inside and outside the focus
+	in_area_density.append(incount * np.mean(insize) / fsarea)
+	out_area_density.append(outcount * np.mean(outsize) / (1. - fsarea))
+	                        
+	# Point densities inside and outside the focus
+	in_point_density.append(incount / fsarea)
+	out_point_density.append(outcount / (1. - fsarea))
+	print "Densities calculated."
+
+
+
+
+
+
+
+
+
+
+
+
+
+	# Variance and CoV of outside the focus
+	mask = np.zeros_like(clean)
+	for c in cent :
+	    mask[int(np.round(c[0])), int(np.round(c[1]))] = 1
+	mask = ndimage.binary_fill_holes(ndimage.maximum_filter(mask, 35))
+
+	areadensity = filter.gaussian_filter(clean, 51)
+	var_area_density.append(np.var(areadensity[mask == 0]))
+	cov_area_density.append(np.std(areadensity[mask == 0] / np.mean(areadensity[mask == 0])))
+
+	pointdensity = np.zeros_like(clean)
+	for n in nuclei :
+	    pointdensity[int(np.round(n.centroid[0])), int(np.round(n.centroid[1]))] = 1
+	pointdensity = filter.gaussian_filter(pointdensity, 51)
+
+	var_point_density.append(np.var(pointdensity[mask == 0]))
+	cov_point_density.append(np.std(pointdensity[mask == 0] / np.mean(pointdensity[mask == 0])))
+	print "Variances calculated. "
+
+
+
+
+
+
+
+
+
+
+	# Pairwise distances between nuclei outside the focal area 
+	distances = distance.squareform(distance.pdist([n.centroid for n in outer_nuclei]))
+
+	for i in range(len(distances)) :
+	    distances[i, i] = np.inf
+	    
+	pairwise.append(np.min(distances, axis=0).mean())
+	print "Pairwise distances calculated."
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+in_count = []
+out_count = []
+in_size = []
+out_size = []
+in_area_density = []
+out_area_density = []
+in_point_density = []
+out_point_density = []
+var_point_density = []
+cov_point_density = []
+var_area_density = []
+cov_area_density = []
+fs_area = []
+pairwise = []
+
+results = {}
+results["in_count"] = in_count
+results["out_count"] = out_count
+results["in_size"] = in_size
+results["out_size"] = out_size
+results["in_area_density"] = in_area_density
+results["out_area_density"] = out_area_density
+results["in_point_density"] = in_point_density
+results["out_point_density"] = out_point_density
+results["name"] = name
+results["ID"] = imageid
+results["var_point_density"] = var_point_density
+results["cov_point_density"] = cov_point_density
+results["var_area_density"] = var_area_density
+results["cov_area_density"] = cov_area_density
+results["fs_area"] = fs_area
+results["pairwise"] = pairwise
+
+F = open("results/thread_%s" % sys.argv[1], "w")
+pickle.dump(results, F)
+F.close()
+print "Thread %s finished." % sys.argv[1]
+