#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Nov  8 11:22:56 2022

@author: george
"""

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.widgets import Slider
from sklearn.neighbors import KDTree
import random
from tqdm import tqdm
import os, glob
from skimage import io
from skimage.transform import resize
from scipy.spatial import distance_matrix
from sklearn.neighbors import KernelDensity

from sklearn.cluster import DBSCAN
from sklearn import metrics
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler

import math
import numpy as np
from scipy import stats
from scipy import spatial
from scipy import histogram2d

%matplotlib qt 

from distutils.version import StrictVersion
import flika
from flika import global_vars as g
from flika.window import Window
from flika.process.file_ import save_file_gui, open_file_gui

from qtpy.QtGui import QColor

flika_version = flika.__version__
if StrictVersion(flika_version) < StrictVersion('0.2.23'):
    from flika.process.BaseProcess import BaseProcess, WindowSelector, SliderLabel, CheckBox
else:
    from flika.utils.BaseProcess import BaseProcess, WindowSelector, SliderLabel, CheckBox

from flika import *
%gui qt

#recording options
pixelSize = 0.108

def loadPointData(trackFile, filterPoints =False, filterType = None, classType=None):   
  
    ######### load FLIKA pyinsight data into DF ############
    pointsDF = pd.read_csv(trackFile)
    pointsDF['frame'] = pointsDF['frame'].astype(int)
    pointsDF['track_number'] = pointsDF['track_number'].astype(int)
    
    if filterPoints==True:
        #choose filter
        if filterType == 'SVM':
            pointsDF = pointsDF[pointsDF['SVMPredictedClass'] == classType]  
                    
    return pointsDF


def loadImgData(tiffFile, transpose=False, rotateflip=True):
    #load image data
    img = io.imread(tiffFile) 
    
    if transpose == True:
        #reshape to x,y,t
        img = np.transpose(img,(1,2,0))
        
    if rotateflip == True:
        #rotate and flip for flika
        img = img[:,:,::-1]
        img = np.rot90(img, k=1, axes=(1,2))
    
    return img


def plotDataOnStack(img, points, pixelSize, crop=True):
    #plot point data on tiff stack
    
    pointWindow = Window(img)
        
    points_byFrame = points[['frame','x','y']]

    #points_byFrame['point_color'] = QColor(g.m.settings['point_color'])
    #points_byFrame['point_size'] = g.m.settings['point_size']
    pointArray = points_byFrame.to_numpy()
    
    pointWindow.scatterPoints = [[] for _ in np.arange(pointWindow.mt)]
    
    
    for pt in pointArray:
        t = int(pt[0])
        if pointWindow.mt == 1:
            t = 0
        pointSize = g.m.settings['point_size']
        pointColor = QColor(g.m.settings['point_color'])
        #position = [pt[1]+(.5* (1/pixelSize)), pt[2]+(.5* (1/pixelSize)), pointColor, pointSize]
        position = [pt[1], pt[2], pointColor, pointSize]    
        pointWindow.scatterPoints[t].append(position)
    pointWindow.updateindex()
    
    
def plotTracks(tracksDF):
    fig1,ax1 = plt.subplots(1, 1, sharex=True, sharey=True)
               
    groups_points = tracksDF.groupby('track_number')
    for name, group in groups_points:
        ax1.plot(group.x, group.y, marker='.', linestyle='-', markersize=3, label=name)
    
    #ax1.set_title(exptName)        
    #ax1.set_xlim(xaxisMin,xaxisMax)
    #ax1.set_ylim(yaxisMin,yaxisMax)
    
    ax1.invert_yaxis()                
    ax1.set(xlabel='PositionX [pixels]')
    ax1.set(ylabel='PositionY [pixels]')

    plt.show()  




def plotDensityHistogram(tracksDF, filteredTracksDF):
    xdat = tracksDF['x'].to_numpy()
    ydat = tracksDF['y'].to_numpy()    
    
    xfiltered = filteredTracksDF['x'].to_numpy()
    yfiltered = filteredTracksDF['y'].to_numpy()
    
    #histogram definition
    xyrange = [[min(xdat),max(xdat)],[min(ydat),max(ydat)]] # data range
    bins = [1000,1000] # number of bins
    thresh = 2  #density threshold
    
    # histogram the data
    hh, locx, locy = histogram2d(xdat, ydat, range=xyrange, bins=bins)
    posx = np.digitize(xdat, locx)
    posy = np.digitize(ydat, locy)
    
    #select points within the histogram
    ind = (posx > 0) & (posx <= bins[0]) & (posy > 0) & (posy <= bins[1])
    hhsub = hh[posx[ind] - 1, posy[ind] - 1] # values of the histogram where the points are
    #xdat1 = xdat[ind][hhsub < thresh] # low density points
    #ydat1 = ydat[ind][hhsub < thresh]
    hh[hh < thresh] = np.nan # fill the areas with low density by NaNs
    hh[hh > 10] = np.nan # fill the areas with low density by NaNs    
    
    plt.imshow(np.flipud(hh.T),cmap='Set1',extent=np.array(xyrange).flatten(), interpolation='none', origin='upper')
    plt.colorbar()   
    plt.plot(xfiltered, yfiltered, '.',color='blue')
    plt.show()




def kde2D(x, y, bandwidth, xbins=256j, ybins=256j, **kwargs): 
    """Build 2D kernel density estimate (KDE)."""

    # create grid of sample locations (default: 100x100)
    xx, yy = np.mgrid[x.min():x.max():xbins, 
                      y.min():y.max():ybins]

    xy_sample = np.vstack([yy.ravel(), xx.ravel()]).T
    xy_train  = np.vstack([y, x]).T

    kde_skl = KernelDensity(bandwidth=bandwidth, **kwargs)
    kde_skl.fit(xy_train)

    # score_samples() returns the log-likelihood of the samples
    z = np.exp(kde_skl.score_samples(xy_sample))
    return xx, yy, np.reshape(z, xx.shape)


def plotKDE(tracksDF,filteredTracksDF):
    x = tracksDF['x'].to_numpy()
    y = tracksDF['y'].to_numpy()

    xfiltered = filteredTracksDF['x'].to_numpy()
    yfiltered = filteredTracksDF['y'].to_numpy()
    
    xx, yy, zz = kde2D(x, y, 5.0)
    
    plt.pcolormesh(xx, yy, zz)
    plt.scatter(xfiltered , yfiltered , s=2, facecolor='white')    
    
    return xx, yy, zz
     

def dbscan(tracksDF):
    X = np.column_stack([tracksDF['x'].to_list(),tracksDF['y'].to_list()])

    db = DBSCAN(eps=1, min_samples=3).fit(X)
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_
    
    # Number of clusters in labels, ignoring noise if present.
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
    n_noise_ = list(labels).count(-1)
    
    print("Estimated number of clusters: %d" % n_clusters_)
    print("Estimated number of noise points: %d" % n_noise_)
    print("Silhouette Coefficient: %0.3f" % metrics.silhouette_score(X, labels))
    
    return labels, X, core_samples_mask, n_clusters_

def plotDBSCAN(tracksDF):
    
    labels, X, core_samples_mask, n_clusters_ = dbscan(tracksDF)
    
    # Black removed and is used for noise instead.
    unique_labels = set(labels)
    colors = [plt.cm.Spectral(each) for each in np.linspace(0, 1, len(unique_labels))]
    for k, col in zip(unique_labels, colors):
        if k == -1:
            # Black used for noise.
            col = [0, 0, 0, 1]
    
        class_member_mask = labels == k
    
        xy = X[class_member_mask & core_samples_mask]
        plt.plot(
            xy[:, 0],
            xy[:, 1],
            "o",
            markerfacecolor=tuple(col),
            markeredgecolor="k",
            markersize=14,
        )
    
        xy = X[class_member_mask & ~core_samples_mask]
        plt.plot(
            xy[:, 0],
            xy[:, 1],
            "o",
            markerfacecolor=tuple(col),
            markeredgecolor="k",
            markersize=6,
        )
    
    plt.title("Estimated number of clusters: %d" % n_clusters_)
    plt.show()
    return labels, X, core_samples_mask, class_member_mask

if __name__ == '__main__':    
    
    #non-BAPTA
    trackFilePath = '/Users/george/Data/10msExposure10s/wt/GB_165_2022_03_01_HTEndothelial_NonBapta_plate1_1_MMStack_Default_bin10_locs_tracksRG_SVMPredicted.csv'
    tiffFilePath = '/Users/george/Data/10msExposure10s/wt/GB_165_2022_03_01_HTEndothelial_NonBapta_plate1_1_MMStack_Default_bin10.tif'
    
    #non-BAPTA FIXED
    #trackFilePath = '/Users/george/Data/10msExposure10s_fixed/wt/GB_195_2022_08_31_HTEndothelial_NonBAPTA_Fixed_plate1_1_MMStack_Default_bin10_locs_tracksRG_SVMPredicted.csv'
    #tiffFilePath = '/Users/george/Data/10msExposure10s_fixed/wt/GB_195_2022_08_31_HTEndothelial_NonBAPTA_Fixed_plate1_1_MMStack_Default_bin10.tif'
    
    tracksDF = loadPointData(trackFilePath, filterPoints = False, filterType = 'SVM', classType=1)
    filteredTracksDF = loadPointData(trackFilePath, filterPoints = True, filterType = 'SVM', classType=3)
    #img= loadImgData(tiffFilePath, transpose=False, rotateflip=True)    
    
    #fa = start_flika()   
    #plotDataOnStack(img, tracksDF, pixelSize)
    #plotTracks(tracksDF)

    #plotDensityHistogram(tracksDF,filteredTracksDF)
    #plotKDE(tracksDF,filteredTracksDF)
    labels, X, core_samples_mask, class_member_mask = plotDBSCAN(tracksDF)

    
    #fa.close()
    
    
    
    