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]