# -*- 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 sklearn.linear_model import BayesianRidge

from sklearn.mixture import GMM

from scipy import spatial, ndimage, signal, stats

import matplotlib

%matplotlib inline

import matplotlib.pyplot as plt

import seaborn

matplotlib.rcParams["figure.figsize"] = (16, 10)

#%matplotlib inline

# WANT

# - Density outside inflammatory zone

# - Area, width, length of focus

# - Area of inflammation / area of vein

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

image = exposure.equalize_adapthist(io.imread(files[4]))

#image = io.imread(files[4])

io.imshow(image)

plt.grid(False)

# <codecell>

#binary = np.logical_or(filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 0], cutoff=0.4, gain=20), lengthscale), \

#                       filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 2], cutoff=0.5, gain=20), lengthscale))

binary = filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 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>

(xdim, ydim, _) = image.shape

xdim /= 10

ydim /= 10

local_density = filter.gaussian_filter(clean, 61)

fig, ax = plt.subplots()

io.imshow(local_density)

#for n, c in enumerate(measure.find_contours(local_density, local_density.max()/2)) :

#    ax.plot(c[:, 1], c[:, 0], linewidth=2)

plt.grid(False)

# <codecell>

ent = filter.gaussian_filter(filter.rank.entropy(local_density, morphology.disk(3)), 75)

ent -= ent.min()

ent /= ent.max()

local_density -= local_density.min()

local_density /= local_density.max()

io.imshow(ent)

plt.grid(False)

# <codecell>

info = ent * (1 + local_density)

io.imshow(info)

plt.grid(False)

# <codecell>

bw = (info) > filter.threshold_otsu(info)

#bw = morphology.dilation(bw, morphology.disk(21))

fig, ax = plt.subplots()

io.imshow(image)

C = measure.find_contours(bw, 0.5)

centroid = []

vals = []

for c in C :

    ax.plot(c[:, 1], c[:, 0], lw=5)

    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))]

ax.plot(cent[:, 1], cent[:, 0], lw=5, c="k", alpha = 0.7)

plt.grid(False)

# <codecell>

plt.scatter(centroid, centroid / np.array(vals))

#ax = plt.gca()

#ax.set_yscale("log")

#ax.set_ylim([1e-9, 1e-5])

print 1. / np.array(vals)

# <codecell>

# DO NOT RUN ME

for f in files :

    image = exposure.equalize_adapthist(io.imread(f))

    binary = filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 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)

    local_density = filter.gaussian_filter(clean, 61)

    ent = filter.gaussian_filter(filter.rank.entropy(local_density, morphology.disk(3)), 101)

    ent -= ent.min()

    ent /= ent.max()

    local_density -= local_density.min()

    local_density /= local_density.max()

    info = ent * (1 + local_density)

    bw = (info) > filter.threshold_otsu(info)

    fig, ax = plt.subplots()

    plt.imshow(image)

    C = measure.find_contours(bw, 0.5)

    centroid = []

    for c in C :

        centroid.append(np.linalg.norm([c[:, 1].mean() - bw.shape[1] / 2, c[:, 0].mean() - bw.shape[0] / 2]))

    cent = C[np.argmin(centroid)]

    ax.scatter(cent[:, 1], cent[:, 0], lw=3)

    plt.grid(False)

    plt.savefig(f.split("Sheep")[1])

    print ("Done " + f.split("Sheep")[1])

# <codecell>

f.split("Sheep")

# <codecell>

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

# <codecell>

io.imshow(filter.gaussian_filter(image, 31))

# <codecell>

# Local density

B = image[:, :, 2]

#B = filter.gaussian_filter(np.abs(B - B.mean()), 31)

#io.imshow(B > 1.1*filter.threshold_otsu(B))

io.imshow(color.rgb2xyz(image)[:, :, 2])

#local_density = filter.gaussian_filter(image[:,:,0], 41)# - filter.gaussian_filter(clean, 201)

#io.imshow(local_density > 1.2 * filter.threshold_otsu(local_density, nbins = 1000))

#io.imshow(local_density - local_density.mean() > 0)

"""

#io.imshow(filter.gaussian_filter(clean, 51))

Q = (signal.convolve2d(ski.img_as_float(clean / clean.max()), ski.img_as_float(morphology.disk(201)), "valid"))# - filter.gaussian_filter(clean, 251))

Q -= Q.min()

Q /= Q.max()

io.imshow(Q)

"""

# <codecell>

B = signal.fftconvolve(ski.img_as_float(A), morphology.disk(31))

B -= B.min()

B /= B.max()

isignal.fftconvolve

# <codecell>

import numpy as np

import matplotlib.pyplot as plt

from scipy import signal

B = np.zeros((100,100))

B[30:70, 30:70] = 1

using_c2d = signal.convolve2d(B, B, mode = "same")

using_fft = signal.fftconvolve(B, B, mode = "same")

plt.subplot(131)

plt.imshow(B)

plt.grid(False)

plt.subplot(132)

plt.imshow(using_c2d)

plt.grid(False)

plt.subplot(133)

plt.imshow(using_fft)

plt.grid(False)

print ((using_fft - using_c2d) < 1e-10).all()

# <codecell>

plt.imshow(using_fft - using_c2d)

# <codecell>

def fftconvolve2d(x, y, mode="full"):

    """

    x and y must be real 2-d numpy arrays.

    mode must be "full" or "valid".

    """

    x_shape = np.array(x.shape)

    y_shape = np.array(y.shape)

    z_shape = x_shape + y_shape - 1

    z = ifft2(fft2(x, z_shape) * fft2(y, z_shape)).real

    if mode != "full":

        # To compute a valid shape, either np.all(x_shape >= y_shape) or

        # np.all(y_shape >= x_shape).

        valid_shape = x_shape - y_shape + 1

        if np.any(valid_shape < 1):

            valid_shape = y_shape - x_shape + 1

            if np.any(valid_shape < 1):

                raise ValueError("empty result for valid shape")

        if mode == "valid" :                    

            start = (z_shape - valid_shape) // 2

            end = start + valid_shape

            z = z[start[0]:end[0], start[1]:end[1]]

        elif mode == "same" :

            start = (z_shape - x_shape) // 2

            z = z[start[0]:-start[0], start[1]:-start[1]]

    return z

def makeGaussian(size, fwhm = 3, center=None):

    """ Make a square gaussian kernel.

    size is the length of a side of the square

    fwhm is full-width-half-maximum, which

    can be thought of as an effective radius.

    """

    x = np.arange(0, size, 1, float)

    y = x[:,np.newaxis]

    if center is None:

        x0 = y0 = size // 2

    else:

        x0 = center[0]

        y0 = center[1]

    return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)

A = clean[1000:2000, 2500:]

io.imshow(fftconvolve2d(A, makeGaussian(61, fwhm=31), mode="same"))

# <codecell>

# <codecell>

for im in files :

    # Read image

    image = io.imread(im)

    # Threshold for nuclei

    binary = np.logical_or(filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 0], cutoff=0.4, gain=30), lengthscale), \

                      filter.threshold_adaptive(exposure.adjust_sigmoid(image[:, :, 2], gain=30), lengthscale) )

    # Clean up binary image

    clean = (1 - morphology.binary_closing(morphology.remove_small_objects(binary, lengthscale), morphology.disk(3))) #* \

#            color.rgb2grey(image)

    io.imshow(image)

    # Create smooth density surface

#    surf = filter.gaussian_filter(clean, lengthscale, mode = "constant")

    # Calculated baseline and remove

    """

    r = [BayesianRidge().fit(np.expand_dims(np.arange(surf.shape[i]), axis=1), np.mean(surf, axis = 1-i)) for i in [0, 1]]

    X, Y = np.meshgrid(r[1].predict(np.expand_dims(np.arange(surf.shape[1]), axis=1)), \

                       r[0].predict(np.expand_dims(np.arange(surf.shape[0]), axis=1)))

    """

"""

    r = [BayesianRidge().fit(np.vander(np.arange(surf.shape[i]), 2), np.mean(surf, axis = 1-i)) for i in [0, 1]]

    X, Y = np.meshgrid(r[1].predict(np.vander(np.arange(surf.shape[1]), 2)), \

                       r[0].predict(np.vander(np.arange(surf.shape[0]), 2)))

    baseline = X * Y

#    hist, bins = np.histogram(surf.ravel(), 100)

#    baseline = bins[np.argmax(hist)+1]# + np.std(surf)

    zeroed = surf - baseline

    zeroed[zeroed < 0] = 0

    # Find peaks in image centre

    peaks = feature.peak_local_max(zeroed[offset:-offset, offset:-offset])

    dist = np.min(spatial.distance.pdist(peaks))/2. if len(peaks) > 1 else min(zeroed.shape)/2.

    # Compute volume integrals and maxima per peak

    integrals = []

    maxima = []

    for i in peaks :

        mask = np.zeros_like(zeroed)

        mask[draw.circle(i[0] + offset, i[1] + offset, dist, shape = zeroed.shape)] = 1

        integrals.append(np.sum(zeroed * mask) / np.sum(mask))

        maxima.append((zeroed[i[0] + offset, i[1] + offset] + baseline) / baseline)

"""    

#    io.imshow(surf)

#    plt.grid(0)

# <codecell>

"""

#image_max = ndimage.maximum_filter(dist, size=11, mode='constant')

#dist = ndimage.distance_transform_edt(clean)

#crit = np.zeros_like(clean)

#coordinates = feature.peak_local_max(dist, min_distance=5)

s2 = filter.gaussian_filter(clean, 21)

dist2 = ndimage.distance_transform_edt(s2 > filter.threshold_otsu(s2))

coordinates2 = feature.peak_local_max(dist2, min_distance=5)

#print len(coordinates), len(coordinates2)

#scales = [11, 21, 31, 41, 51]

spatial.voronoi_plot_2d(spatial.Voronoi(coordinates2))

#io.imshow(filter.gaussian_filter(clean, 11))

"""

nuclei = feature.peak_local_max(ndimage.distance_transform_edt(clean), min_distance=11)

tessellation 

scales = [11, 21, 31, 41, 51]

for i in scales :

# <codecell>

x = []

for i in range(3, 51, 2) :

    s = filter.gaussian_filter(clean, i)

    x.append(len(feature.peak_local_max(ndimage.distance_transform_edt(s > filter.threshold_otsu(s)), min_distance=5)))

plt.plot(range(3, 51, 2), x)

# <codecell>

Q = [clean]

kernelsize = 51

downsamplesize = 3

for i in range(1, 6) :

    im = filter.gaussian_filter(morphology.binary_opening(Q[i-1], morphology.disk(5)), kernelsize)

    im2 = zeros((im.shape[0]/downsamplesize, im.shape[1]/downsamplesize))

    for x in range(im2.shape[0]) :

        for y in range(im2.shape[1]) :

            im2[x, y] = np.mean(im[downsamplesize*x:downsamplesize*(x+1), downsamplesize*y:downsamplesize*(y+1)])

    im2 -= im2.min()

    im2 /= im2.max()

    Q.append(np.round(im2*1.))

# <codecell>

q = []

for i in Q :

    q.append(np.sum(i) / (i.shape[0] * i.shape[1]))

plt.plot(q)

# <codecell>

io.imshow(Q[2])

# <codecell>

#io.imshow(morphology.binary_erosion(clean, np.ones((21,21))))

#io.imshow(morphology.binary_dilation(clean, morphology.disk(25)))

io.imshow(morphology.binary_dilation(filter.rank.mean(clean, morphology.disk(51)), morphology.disk(51)))

# <codecell>

blah = filter.gaussian_filter(clean, 401)

plt.hist(blah.ravel(), 256);

# <codecell>

io.imshow(blah)# > np.median(blah.ravel()))

# <codecell>

#r = [BayesianRidge().fit(np.expand_dims(np.arange(surf.shape[i]), axis=1), np.mean(surf, axis = 1-i)) for i in [0, 1]]

X, Y = np.meshgrid(r[0].predict(np.expand_dims(np.arange(surf.shape[0]), axis=1)), \

                   r[1].predict(np.expand_dims(np.arange(surf.shape[1]), axis=1)))

zeroed.shape

# <codecell>

c = clean - np.min(clean)

c /= c.max()

c = c.astype(bool)

io.imsave("/Users/qcaudron/Desktop/charo/2_smoothed.jpg", ski.img_as_uint(surf))

# <codecell>

z1 = np.mean(surf, axis=0)

z2 = np.mean(surf, axis=1)

#for i in range(surf.shape[1]) : 

#    plt.plot(surf[:, i], "k")

#plt.plot(z2)

r = [BayesianRidge().fit(np.vander(np.arange(surf.shape[i]), 2), np.mean(surf, axis = 1-i)) for i in [0, 1]]

r1 = BayesianRidge().fit(np.arange(len(z1)).reshape(len(z1),1), z1)

r2 = BayesianRidge().fit(np.arange(len(z2[500:-500])).reshape(len(z2[500:-500]),1), z2[500:-500])

#plt.plot(r1.predict(np.arange(len(z1)).reshape(len(z1),1)), linewidth=5)

plt.plot(r2.predict(np.arange(len(z2)).reshape(len(z2),1)), linewidth=5)

plt.plot(z2, linewidth=5)

#plt.axhline(b[np.argmax(h)], c="r", linewidth=3)

#plt.plot(r[0].predict(np.vander(np.arange(surf.shape[0]), 2)), linewidth=3)

#plt.plot(r[0].predict(np.arange(len(z1)).reshape(len(z1),1)), linewidth=3)

#plt.plot(r[0].predict(np.expand_dims(np.arange(surf.shape[0]), axis=1)), linewidth=5)

#plt.axhline(np.mean(z1 / r1.predict(np.arange(len(z1)).reshape(len(z1),1))))

# <codecell>

lz = np.log(z2)

r3 = BayesianRidge().fit(np.arange(len(lz[500:-500])).reshape(len(lz[500:-500]),1), lz[500:-500])

plt.plot(np.exp(lz))

plt.plot(np.exp(r3.predict(np.arange(len(lz)).reshape(len(lz),1))))

# <codecell>

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

x, y = np.meshgrid(range(surf.shape[0]), range(surf.shape[1]))

x.shape

ax.plot_surface(y, x, zeroed.T, rstride = 50, cstride = 50, antialiased = False, cmap = cm.coolwarm, linewidth=0)

# <codecell>

np.min(baseline)

# <codecell>

x2, y2 = draw.circle(peaks[1][0]+1000, peaks[1][1]+1000, dist*1., shape = (zeroed.shape))

mask2 = zeros_like(zeroed)

mask2[x, y] = 1

np.sum(zeroed * mask - zeroed * mask2)

np.sum(zeroed * mask2)