#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