# -*- 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