#Warnings filter
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning) 
warnings.filterwarnings("ignore", category=RuntimeWarning) 

#General library
import numpy as np
import matplotlib.pyplot as plt

#Path
import os, os.path
import pickle
import shutil
import glob

# Interactive input
from tkinter.filedialog import askopenfilename
from tkinter.filedialog import askdirectory

#Modeling
from sklearn.cluster import KMeans 
from sklearn.cluster import DBSCAN 
from sklearn.decomposition import PCA
from sklearn import svm, datasets
from sklearn.metrics import auc
from sklearn.metrics import plot_roc_curve
from sklearn.model_selection import StratifiedKFold
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report
from sklearn.metrics import plot_confusion_matrix
from sklearn.metrics import confusion_matrix

# Scipy
from scipy.misc import face
from scipy import ndimage as ndi
from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
from skimage.segmentation import clear_border, mark_boundaries
from skimage.filters import roberts, sobel
from scipy.signal.signaltools import wiener

#Image
import cv2
from PIL import Image
import matplotlib.pyplot as plt

#Medical Image
import pydicom
from pydicom import dcmread
from pydicom.data import get_testdata_file
from roipoly import RoiPoly
import mahotas
import mahotas.demos
import mahotas as mh
import numpy as np
from pylab import imshow, show

# ipywidgets for some interactive plots
from ipywidgets.widgets import * 
import ipywidgets as widgets

# Load the model
pickle_in = open('svc_pca_classifier.pickle',"rb")
model = pickle.load(pickle_in)

# Load the PCA object
pickle_in = open('pca.pickle',"rb")
pca = pickle.load(pickle_in)

# Function for read and selects three scans from a folder exam
def select_paths(path):
    lista = glob.glob(path + "\*.dcm")
    list_idx = []
    n_dcm = len(lista)
    if n_dcm > 50:
        n_dcm_interval = int(n_dcm / 12)
    else:
        print("The selected folder has less than 50 images.. try to select another folder that contain more scans")
    lista_idx = [n_dcm_interval*4, n_dcm_interval*5, n_dcm_interval*6]
    return lista, lista_idx

# This function extract and print some metadata related to the patient
def get_metadata(folder_path):
    lista = glob.glob(folder_path + "\\" + "*.dcm")
    n = len(lista)
    ds = dcmread(lista[0])
    print("PatientID: " + str(ds.PatientID))
    print("\tSex: " + str(ds.PatientSex))
    print("\tAge: " + str(int(ds.PatientAge[:-1])))
    print("\nStudy Description: " + str(ds.StudyDescription))
    print("\tExam Modality: " + str(ds.Modality))
    print("\tNumber of scans: " + str(len(lista)))
    print("\tImage size " + str(ds.Rows) + "x" + str(ds.Columns))
    print("\tPatient Position: " + str(ds.PatientPosition))
    print("\tSlice Thickness: " + str(ds.SliceThickness))
    print("\tSpace between slices: " + str(ds.SpacingBetweenSlices))
    return None

# This function load and process the image related to the specified path
def load_and_preprocess_dcm(path):
    ds = dcmread(path)
    img = ds.pixel_array
    img = ds.RescaleIntercept + img * ds.RescaleSlope
    minimum = np.amin(np.amin(img))
    maximum = np.amax(np.amax(img))
    img[img == minimum] = -1000
    # Normalization
    img = (img - minimum) / maximum
    new_img = ((img - img.min()) * (1/(img.max() - img.min()) * 255)).astype('uint8')
    new_img = new_img[50:460,30:480]
    return new_img

# A simple function to show the three images selected in a subplot
def select_and_show_3(lista, lista_idx):
    fig, ([ax1, ax2, ax3]) = plt.subplots(1,3, figsize=(15,12))
    img1 = load_and_preprocess_dcm(lista[lista_idx[0]])
    img2 = load_and_preprocess_dcm(lista[lista_idx[1]])
    img3 = load_and_preprocess_dcm(lista[lista_idx[2]])
    ax1.imshow(img1, cmap=plt.cm.gray)
    ax1.set_title("Scansione n: " + str(lista_idx[0]))
    ax2.imshow(img2, cmap=plt.cm.gray)
    ax2.set_title("Scansione n: " + str(lista_idx[1]))
    ax3.imshow(img3, cmap=plt.cm.gray)
    ax3.set_title("Scansione n: " + str(lista_idx[2]))
    return [img1, img2, img3]

# A simple function to show the three mask in a subplot
def show_3_mask(mask_list, idx_sel):
    fig, ([ax1, ax2, ax3]) = plt.subplots(1,3, figsize=(15,12))
    ax1.imshow(mask_list[0], cmap=plt.cm.gray)
    ax1.set_title("Maschera della scansione n: " + str(idx_sel[0]))
    ax2.imshow(mask_list[1], cmap=plt.cm.gray)
    ax2.set_title("Maschera della scansione n: " + str(idx_sel[1]))
    ax3.imshow(mask_list[2], cmap=plt.cm.gray)
    ax3.set_title("Maschera della scansione n: " + str(idx_sel[2]))
    return None

# This function allow to create a lung mask with the kmeans clustering algorithm
def mask_kmeans(img_array):
    img_reshaped = img_array.reshape(-1,1)
    kmeans = KMeans(n_clusters = 2, random_state = 99).fit(img_reshaped)
    kmeans_y = kmeans.fit_predict(img_reshaped)
    lab = kmeans.labels_
    mask = lab.reshape(img_array.shape)
    # Control the mask
    if mask[1,1] == 0:
        mask = abs(mask-1)
    return mask

# Function to make the mask more robust
def fill_mask(first_mask):
    cleared = clear_border(first_mask)
    selem = disk(1.5)
    binary = binary_erosion(cleared, selem)
    selem = disk(10)
    binary = binary_closing(binary, selem)
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
    final_mask = binary.astype('uint8')
    # Control the mask
    n_0 = len(final_mask[final_mask==0])
    n_1 = len(final_mask[final_mask==1])
    if (n_1/(n_0 + n_1)) < 0.1 or (n_1/(n_0 + n_1)) > 0.9:
        print("The selected scan is not good enough")
    return final_mask

# Function to manually select the lungh mask or any other ROI
def manually_roi_select(lista_img, mask_list, n_mask):
    fig = plt.figure()
    plt.imshow(lista_img[int(n_mask)-1], interpolation='nearest', cmap=plt.cm.gray)
    plt.colorbar()
    plt.title("left click: line segment right click or double click: close region")
    roi1 = RoiPoly(color='r', fig=fig)
    mask = roi1.get_mask(mask_list[int(n_mask)-1])
    mask_list[int(n_mask)-1] = mask
    return roi1, mask

# Function to show the manually selected ROI
def show_roi_selected(roi1, mask, n_mask, lista_img):
    fig, ([ax1, ax2]) = plt.subplots(1,2, figsize=(15,12))
    roi1.display_roi()
    ax1.imshow(mask, cmap=plt.cm.gray)
    ax1.set_title("Maschera")
    ax2.imshow(lista_img[int(n_mask)-1], cmap=plt.cm.gray)
    ax2.set_title("Maschera sull'immagine originale")
    return None

# This function take the image and extract the Haralick Feature matrix (4x13)
def extract_haralick(img):
    # adding gaussian filter
    img_test = mahotas.gaussian_filter(img, .5, cval = 100, mode = "reflect")
    # setting threshold (threshed is the mask)
    threshed = (img_test > np.quantile(img_test[img_test!=0],.95))
    # making is labeled image
    feature_map, n = mahotas.label(threshed)
    # getting haralick features
    h_feature = mahotas.features.haralick(feature_map, distance = 2)    
    return feature_map, h_feature

# Function for extract and show the equalized and filtered image (dst) and the Haralick map
# for the three scans selected for a patient
def features_and_plot(lista_img, lista_mask, printing=False):
    
    id_feat_list = []
    row, col = lista_img[0].shape
    
    for cont in range(0,len(lista_img)):
        # Equalized hist of image
        dst = cv2.equalizeHist(lista_img[cont]*lista_mask[cont])
        # Wiener filter
        filtered_img = wiener(dst, (3, 3), noise=10)
        # Extract Haralick Features
        feature_map, h_feature = extract_haralick(filtered_img)
        id_feat_list.append(h_feature)

        # showing the features
        if printing:
            fig, ([ax1, ax2]) = plt.subplots(1,2, figsize=(12,10))
            fig.suptitle("Scans: " + str(cont + 1), y = 0.8, fontsize=16)
            ax1.imshow(dst, cmap=plt.cm.gray)
            ax1.set_title("Equalized image (Masked)")
            ax2.imshow(feature_map, cmap=plt.cm.viridis)
            ax2.set_title("Haralick Features Map")
            plt.tight_layout()         
    return id_feat_list   

# Function to predict the probability that the patient is positive or negative to Covid-19
def covid_predict(feat):
    lista_pat=[]
    for elem in feat:
        lista_pat.append(elem.reshape(-1))
    pat_array = np.vstack(lista_pat)
    pat_array = np.reshape(lista_pat, (156))
    pat_array = pat_array.reshape(1,-1)
    input_model = pca.transform(pat_array)
    idx_predicted = model.predict(input_model)
    proba = model.predict_proba(input_model)[0,idx_predicted]
    proba = np.around(proba,4) * 100
    if idx_predicted == 0:
        print("The patient is negative to Covid-19 with a confidence of " + str(proba[0]) + "%")
    else:
        print("The patient is positive to Covid-19 with a confidence of " + str(proba[0]) + "%")
    return None