#!/usr/bin/env python2
# -*- coding: utf-8 -*-
import tifffile as Tiff
import numpy as np
import matplotlib.pylab as plt
from scipy.signal import medfilt

from skimage.filters import threshold_adaptive, gaussian
from skimage.morphology import closing, disk
from skimage.morphology import remove_small_holes, remove_small_objects
from skimage.feature import canny, blob_log
from skimage.measure import label
from skimage import img_as_float

# Reading the image stacks
# stack1 = img_as_float(Tiff.imread('path_to_image_stack_1.tif'))
# stack2 = ...

# Reading the image stacks
# membranes = []
# ROI1 = (x, y, w, h)
# membranes.append( (stack1, ROI1) )
# ROI2 = ...
# membranes.append( ... )

DDATA = []
for stack, recs in membranes:
    
    # Choose different ROIs
    rec_x, rec_y, rec_w, rec_h = recs
    membrane = stack[:, rec_y:rec_y+rec_h, rec_x:rec_x+rec_w].copy()
    
    # Get dimensions and initialize empty array
    Nslices, H, W = membrane.shape
    Y, X = np.indices((H, W))
    classified = np.zeros_like(membrane, dtype=np.double)
    
    for idx, sslice in enumerate(membrane):
            
            # Remove background
            bkg = gaussian(sslice, 8, mode='mirror')
            frame = sslice - bkg
            frame += np.abs(frame.min())
            
            # Threshold image locally
            threshed = threshold_adaptive(frame, block_size=201)
            threshed = remove_small_holes(threshed, 5)
            
            # Find boundaries and combine both binarized images
            boundary = np.bitwise_xor(threshed, canny(frame, sigma=1.25))
            
            # Morphology operations
            proccd = closing(boundary, disk(1))
            proccd = closing(proccd, disk(1))
            proccd = remove_small_holes(proccd, 30)
            proccd = remove_small_objects(proccd, 10)
            
            # Find small objects in inverted image ...
            blobs = blob_log(1.0-frame, 1.5, 3, threshold=0.05)
            for y, x, s in blobs:
                if proccd[int(y), int(x)] == False:
                    continue
                else:
                    # ... and combine with binarized image
                    iid = (Y-y)**2+(X-x)**2 < (np.sqrt(2)*s)**2
                    proccd[iid] = False
            
            # Store results in empty array
            classified[idx] = np.double(proccd).copy()                    
        
    DDATA.append((classified, membrane))

# To check the segmentation algorithm interactively use the following lines
# in "%pylab qt" mode (in ipython console)
# 
# from skimage.segmentation import find_boundaries
# fig, ax, img = Tiff.imshow(membrane, cmap=plt.cm.rainbow)
# bound = np.asarray([find_boundaries(f) for f in classified])
# bound[bound==True] = np.nan
# Tiff.imshow(bound, figure=fig, alpha=0.2, cmap=plt.cm.gray)

# Physical quantities
um_per_pixel = 0.1582 #µm
dt = 5.0/60.0 #min 

# Visualization
plt.figure(figsize=(6,3))
plt.subplot(122)
for uui, (_bin, _mem) in enumerate(DDATA):
    
    time = np.arange(len(_bin))*dt
    h, w = _bin.shape[1:]
    
    nrSpots = [label(i, background=1).max() for i in _bin]
         
    plt.plot(time, medfilt(np.double(nrSpots)/nrSpots[0]),'-', label=uui)

plt.hlines(1, time.min(), time.max(), linestyles='dashed')
plt.xlabel('time [min]')
plt.ylabel('Normalized domain number')

plt.subplot(121)
for _bin, _mem in DDATA:
    
    time = np.arange(len(_bin))*dt
    h, w = _bin.shape[1:]
    
    ld = _bin.sum(axis=(-2,-1))
    lo = (h*w)-ld
         
    plt.plot(time, medfilt((lo-lo[0])/lo[0]),'-')

plt.xlabel('time [min]')
plt.ylabel(u'$\Delta(A_{Lo})/A_{Lo}$')
plt.ylim([-0.2, 0.2])
plt.yticks([-.2, -.1, 0., .1, .2])

plt.tight_layout()
plt.show()

