{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Part 2 - Tradition Image processing for classification: model training\n\nContinuation of my [previous notebook](https://www.kaggle.com/suryathiru/1-tradition-image-processing-feature-extraction/) where features were manually extracted from the dataset","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:23:59.031199Z","start_time":"2019-11-05T10:23:42.138877Z"},"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport cv2\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom sklearn.model_selection import train_test_split\nimport os\nimport pydicom as dcm","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:23:59.057167Z","start_time":"2019-11-05T10:23:59.043431Z"},"trusted":true},"cell_type":"code","source":"np.random.seed(42)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prepare dataset","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:23:59.426980Z","start_time":"2019-11-05T10:23:59.066494Z"},"trusted":true},"cell_type":"code","source":"data = pd.read_csv('../input/1-tradition-image-processing-feature-extraction/img_features.csv')\ndata.head()","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:23:59.931897Z","start_time":"2019-11-05T10:23:59.431837Z"},"trusted":true},"cell_type":"code","source":"feats = data.features.apply(lambda x: list(eval(x)))  # oops!\ndataset = pd.DataFrame(feats.values.tolist(), \n                        columns=['mean', 'stddev', 'area', 'perimeter', 'irregularity', 'equiv_diam', 'hu1', 'hu2', 'hu4', 'hu5', 'hu6'],\n                       index=data.index)\ndataset['label'] = data['target']\ndataset.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train test split","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:00.079067Z","start_time":"2019-11-05T10:23:59.935607Z"},"trusted":true},"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(dataset.drop('label', axis=1), dataset['label'],\n                                                   test_size=0.4, stratify=dataset['label'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train models\n\nWe evaluate on the following metrics\n1. Accuracy\n2. Precision\n3. Recall\n4. ROC AUC","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:00.098278Z","start_time":"2019-11-05T10:24:00.082231Z"},"trusted":true},"cell_type":"code","source":"from sklearn.metrics import accuracy_score, roc_auc_score, precision_score, recall_score, f1_score\n\ndef print_metrics(y_pred, y_train, yt_pred, y_test):\n    print('Train data metrics:')\n    print('Accuracy: ', accuracy_score(y_train, y_pred))\n    print('Precison: ', precision_score(y_train, y_pred))\n    print('Recall: ', recall_score(y_train, y_pred))\n    print('F1 score: ', f1_score(y_train, y_pred))\n    print('ROC AUC score: ', roc_auc_score(y_train, y_pred))\n    print()\n    print('Test data metrics:')\n    print('Accuracy: ', accuracy_score(y_test, yt_pred))\n    print('Precison: ', precision_score(y_test, yt_pred))\n    print('Recall: ', recall_score(y_test, yt_pred))\n    print('F1 score: ', f1_score(y_test, yt_pred))\n    print('ROC AUC score: ', roc_auc_score(y_test, yt_pred))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Logistic regression model","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:00.979671Z","start_time":"2019-11-05T10:24:00.102522Z"},"trusted":true},"cell_type":"code","source":"from sklearn.linear_model import LogisticRegression\n\nlr = LogisticRegression()\nlr.fit(X_train, y_train)\n\nprint_metrics(lr.predict(X_train), y_train, lr.predict(X_test), y_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Random forest","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:06.306800Z","start_time":"2019-11-05T10:24:00.984193Z"},"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\nrf = RandomForestClassifier(500, max_depth=8, min_samples_split=2,\n                            n_jobs=-1)\nrf.fit(X_train, y_train)\n\nprint_metrics(rf.predict(X_train), y_train, rf.predict(X_test), y_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Gradient boosting classifier","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:23.227137Z","start_time":"2019-11-05T10:24:06.314951Z"},"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import GradientBoostingClassifier\ngb = GradientBoostingClassifier(n_estimators=500, learning_rate=0.01,\n                               max_depth=7, min_samples_split=5)\ngb.fit(X_train, y_train)\n\nprint_metrics(gb.predict(X_train), y_train, gb.predict(X_test), y_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Support Vector Machines","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:30.917546Z","start_time":"2019-11-05T10:24:23.229296Z"},"trusted":true},"cell_type":"code","source":"from sklearn.svm import SVC\nsvm = SVC()\nsvm.fit(X_train, y_train)\n\nprint_metrics(svm.predict(X_train), y_train, svm.predict(X_test), y_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### KNN","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:31.661175Z","start_time":"2019-11-05T10:24:30.920949Z"},"trusted":true},"cell_type":"code","source":"from sklearn.neighbors import KNeighborsClassifier\nknn = KNeighborsClassifier(20)\nknn.fit(X_train, y_train)\n\nprint_metrics(knn.predict(X_train), y_train, knn.predict(X_test), y_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualize the best model","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:33.109109Z","start_time":"2019-11-05T10:24:31.663859Z"},"trusted":true},"cell_type":"code","source":"from yellowbrick.classifier import ConfusionMatrix\n\ncm = ConfusionMatrix(gb, classes=[0,1], \n                     encoder={0: 'normal', 1: 'pneumonia'})\ncm.fit(X_train, y_train)\ncm.score(X_test, y_test)\ncm.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### ROC","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:33.581399Z","start_time":"2019-11-05T10:24:33.112096Z"},"trusted":true},"cell_type":"code","source":"from yellowbrick.classifier import ROCAUC\n\nvisualizer = ROCAUC(gb, classes=[\"normal\", \"pneumonia\"])\n\nvisualizer.fit(X_train, y_train) \nvisualizer.score(X_test, y_test) \nvisualizer.show()    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Feature importance","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:33.952273Z","start_time":"2019-11-05T10:24:33.584565Z"},"trusted":true},"cell_type":"code","source":"importances = rf.feature_importances_\nstd = np.std([tree.feature_importances_ for tree in rf.estimators_],\n             axis=0)\nindices = np.argsort(importances)[::-1]\n\nprint(\"Feature ranking:\")\n\nfor f in range(X_train.shape[1]):\n    print(\"%d. feature %d (%f)\" % (f + 1, indices[f], importances[indices[f]]))\n\nplt.figure(figsize=(13,7))\nplt.title(\"Feature importances\")\n\nplt.bar(range(X_train.shape[1]), importances[indices],\n       color=\"r\", yerr=std[indices], align=\"center\")\nplt.xticks(range(X_train.shape[1]), X_train.columns[indices])\nplt.xlim([-1, X_train.shape[1]])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict with an image","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-10-30T18:37:10.317218Z","start_time":"2019-10-30T18:37:10.312980Z"}},"cell_type":"markdown","source":"### Functions","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:37.368851Z","start_time":"2019-11-05T10:24:33.954992Z"},"trusted":true},"cell_type":"code","source":"import skimage\nfrom skimage import feature, filters","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:37.396732Z","start_time":"2019-11-05T10:24:37.372998Z"},"trusted":true},"cell_type":"code","source":"PATH = '../input/rsna-pneumonia-detection-challenge'\n\ndef load_image(path):    \n    patientImage = path + '.dcm'\n    imagePath = os.path.join(PATH,\"stage_2_train_images/\", patientImage)\n    img = dcm.read_file(imagePath).pixel_array\n    return img\n\ndef imshow_gray(img):\n    plt.figure(figsize=(12,7))\n    return plt.imshow(img, cmap='gray')\n\ndef area(img):\n    # binarized image as input\n    return np.count_nonzero(img)\n\ndef perimeter(img):\n    # edges of the image as input\n    return np.count_nonzero(img)\n\ndef irregularity(area, perimeter):\n    # area and perimeter of the image as input, also called compactness\n    I = (4 * np.pi * area) / (perimeter ** 2)\n    return I\n\ndef equiv_diam(area):\n    # area of image as input\n    ed = np.sqrt((4 * area) / np.pi)\n    return ed\n\ndef get_hu_moments(contour):\n    # hu moments except 3rd and 7th (5 values)\n    M = cv2.moments(contour)\n    hu = cv2.HuMoments(M).ravel().tolist()\n    del hu[2]\n    del hu[-1]\n    log_hu = [-np.sign(a)*np.log10(np.abs(a)) for a in hu]\n    return log_hu\n\n\ndef extract_features(img):\n    mean = img.mean()\n    std_dev = img.std()\n    \n    # hist equalization\n    equalized = cv2.equalizeHist(img)\n    \n    # sharpening\n    hpf_kernel = np.full((3, 3), -1)\n    hpf_kernel[1,1] = 9\n    sharpened = cv2.filter2D(equalized, -1, hpf_kernel)\n    \n    # thresholding\n    ret, binarized = cv2.threshold(cv2.GaussianBlur(sharpened,(7,7),0),0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n    \n    # edge detection\n    edges = skimage.filters.sobel(binarized)\n    \n    # moments from contours\n    contours, hier = cv2.findContours((edges * 255).astype('uint8'),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)\n    select_contour = sorted(contours, key=lambda x: x.shape[0], reverse=True)[0]\n    \n    \n    # feature extraction\n    ar = area(binarized)\n    per = perimeter(edges)\n    irreg = irregularity(ar, per)\n    eq_diam = equiv_diam(ar)\n    hu = get_hu_moments(select_contour)\n    \n    return (mean, std_dev, ar, per, irreg, eq_diam, *hu)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load images","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:43:31.797215Z","start_time":"2019-11-05T10:43:31.461131Z"},"scrolled":false,"trusted":true},"cell_type":"code","source":"test_img = data[data['target']==1]['patientId'].sample(1)\n\nimg = load_image(test_img.values[0])\nimshow_gray(img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Make prediction","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:53.791314Z","start_time":"2019-11-05T10:24:53.715157Z"},"scrolled":true,"trusted":true},"cell_type":"code","source":"feats = list(extract_features(img))\n\npred = gb.predict([feats])\n\nif pred == 1:\n    print('Patient is infected with pneumonia')\nelse:\n    print('Patient is normal')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test on a noisy image","execution_count":null},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:54.300679Z","start_time":"2019-11-05T10:24:53.795055Z"},"trusted":true},"cell_type":"code","source":"from skimage.util import random_noise\n\nimg = random_noise(img, mode='gaussian')\nimg = (img*255).astype('int')\nimshow_gray(img)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-11-05T10:24:54.396406Z","start_time":"2019-11-05T10:24:54.302519Z"},"trusted":true},"cell_type":"code","source":"feats = list(extract_features(img.astype('uint8')))\n\npred = gb.predict([feats])\n\nif pred == 1:\n    print('Patient is infected with pneumonia')\nelse:\n    print('Patient is normal')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python ML","language":"python","name":"ml"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.3"},"varInspector":{"cols":{"lenName":16,"lenType":16,"lenVar":40},"kernels_config":{"python":{"delete_cmd_postfix":"","delete_cmd_prefix":"del ","library":"var_list.py","varRefreshCmd":"print(var_dic_list())"},"r":{"delete_cmd_postfix":") ","delete_cmd_prefix":"rm(","library":"var_list.r","varRefreshCmd":"cat(var_dic_list()) "}},"types_to_exclude":["module","function","builtin_function_or_method","instance","_Feature"],"window_display":false}},"nbformat":4,"nbformat_minor":4}