--- a
+++ b/10x/reference/batch.py
@@ -0,0 +1,206 @@
+# -*- coding: utf-8 -*-
+# <nbformat>3.0</nbformat>
+
+# <codecell>
+
+import os
+
+import numpy as np
+
+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
+
+import matplotlib
+%matplotlib inline
+import matplotlib.pyplot as plt
+import seaborn
+matplotlib.rcParams["figure.figsize"] = (16, 10)
+
+# <codecell>
+
+# Load filenames
+dirq = "/Users/qcaudron/repositories/Quantitative-Histology/10x/data/"
+files = []
+for i in os.listdir(dirq) :
+    if i.endswith(".jpg") :
+        files.append(dirq + i)
+        
+
+offset = 1000
+lengthscale = 301
+A = io.imread(files[4])
+
+
+# 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 = []
+
+
+image = exposure.equalize_adapthist(A)
+io.imshow(image)
+plt.grid(False)
+
+# <codecell>
+
+# Process for nuclei
+
+binary = filter.threshold_adaptive(exposure.adjust_sigmoid(A[:, :, 0], cutoff=0.4, gain = 30), 301).astype(bool)
+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)
+
+io.imshow(clean)
+plt.grid(False)
+
+# <codecell>
+
+# 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))]
+path = matplotlib.path.Path(cent)
+
+io.imshow(image)
+plt.plot(cent[:, 1], cent[:, 0], lw=5, c="k", alpha = 0.7)
+plt.grid(False)
+
+# <codecell>
+
+# Segmentation
+
+distance = ndimage.distance_transform_edt(clean)
+peaks = feature.peak_local_max(distance, indices=False, labels = clean)
+markers = ndimage.label(peaks)[0]
+labels = morphology.watershed(-distance, markers, mask=clean)
+
+io.imshow(labels, interpolation = "nearest", cmap = plt.cm.spectral)
+plt.grid(False)
+
+# <codecell>
+
+# 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))
+
+
+
+
+# Using the shoelace algorithm, compute the area of the focus
+x = path.vertices[:, 0]
+y = path.vertices[:, 1]
+x = np.append(x, path.vertices[0, 0])
+y = np.append(y, path.vertices[0, 1])
+fsarea = np.abs(np.sum(x[:-1] * y[1:] - x[1:] * y[:-1])) / (2. * 3456 * 5184)
+fs_area.append(fsarea)
+
+# <codecell>
+
+# 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))
+
+# <codecell>
+
+# 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])))
+
+# <codecell>
+
+# 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())
+
+# <codecell>
+
+print in_count
+print out_count
+print in_size
+print out_size
+print in_area_density
+print out_area_density
+print in_point_density
+print out_point_density
+print var_point_density
+print cov_point_density
+print var_area_density
+print cov_area_density
+print fs_area
+