{"cells":[{"metadata":{"_uuid":"76f25b25d03405d7ae44958368e878e694942097"},"cell_type":"markdown","source":"# Overview\nHere we take a simple model trained on the miniture tile-based data and apply it to the real competition. \n\nWe follow the following recipe to create the results\n1. Create tiles out of each image\n1. Make a prediction for each tile\n1. Add the square to the results if the prediction score is above a threshold"},{"metadata":{"trusted":true,"_uuid":"4d3df17df66bc4bc1837a00a1c24d5dc39393889"},"cell_type":"code","source":"%matplotlib inline\nimport os\nimport pandas as pd, numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage.io import imread\nimport seaborn as sns # nice visuals\nfrom sklearn.model_selection import train_test_split # splitting data\nfrom matplotlib.patches import Rectangle\nimport pydicom\n# quantifying models\nfrom sklearn.metrics import accuracy_score, classification_report, roc_curve, roc_auc_score, confusion_matrix\ntraining_scale_factor = 2 # how much did we scale the images for training\nDCM_TAG_LIST = ['PatientAge', 'ViewPosition', 'PatientSex']\ndata_dir = '../input/pneumonia-texture-analysis/'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"def categories_to_indicators(in_df):\n    new_df = in_df.copy()\n    new_df['IsMale'] = in_df['PatientSex'].map(lambda x: 'M' in x).astype(float)\n    new_df['IsAP'] = in_df['ViewPosition'].map(lambda x: 'AP' in x).astype(float)\n    return new_df.drop(['PatientSex', 'ViewPosition'], axis=1)\nfull_train_df = categories_to_indicators(pd.read_csv(os.path.join(data_dir, 'train_all.csv')))\nfull_train_stack = imread(os.path.join(data_dir, 'train.tif'))\nfull_train_df.sample(5)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import RobustScaler\nfrom sklearn.pipeline import make_pipeline\n\ndef fit_and_score(in_model, full_features, full_labels, rescale=True):\n    \"\"\"\n    Take a given model, set of features, and labels\n    Break the dataset into training and validation\n    Fit the model\n    Show how well the model worked\n    \"\"\"\n    train_feat, valid_feat, train_lab, valid_lab = train_test_split(full_features, \n                                                                    full_labels,\n                                                                    test_size=0.25,\n                                                                    random_state=2018)\n    \n    if rescale:\n        feature_scaler = RobustScaler()\n        train_feat = feature_scaler.fit_transform(train_feat)\n        valid_feat = feature_scaler.transform(valid_feat)\n    in_model.fit(train_feat, train_lab)\n    predictions = in_model.predict_proba(valid_feat)[:, 1]\n    predicted_class = predictions>0.5\n    tpr, fpr, _ = roc_curve(valid_lab, predictions)\n    auc = roc_auc_score(valid_lab, predictions)\n    print(classification_report(valid_lab, predicted_class))\n    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n    ax1.plot(tpr, fpr, 'r.-', label='Prediction (AUC:{:2.2f})'.format(auc))\n    ax1.plot(tpr, tpr, 'k-', label='Random Guessing')\n    ax1.legend()\n    ax1.set_title('ROC Curve')\n    sns.heatmap(confusion_matrix(valid_lab, predicted_class), \n                annot=True, fmt='4d', ax=ax2)\n    ax2.set_xlabel('Prediction')\n    ax2.set_ylabel('Actual Value')\n    ax2.set_title('Confusion Matrix ({:.1%})'.format(accuracy_score(valid_lab, predicted_class)))\n    return make_pipeline(feature_scaler, in_model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d74e9ea9ec52be9bf1720507927084c403f602d2"},"cell_type":"code","source":"# add the intensity channels\nfull_train_df['Mean_Intensity'] = np.mean(full_train_stack, (1, 2))\nfull_train_df['Std_Intensity'] = np.std(full_train_stack, (1, 2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"47cccaaf1a3948af3a05c3e30e3679c92e2c6378"},"cell_type":"code","source":"from sklearn.svm import SVC\nsvm_model = SVC(probability=True)\nmodel_cols = ['PatientAge', 'opacity_prior', 'IsMale', 'IsAP', \n                   'Mean_Intensity', 'Std_Intensity']\nfitted_pipeline = fit_and_score(\n    svm_model,\n    full_train_df[model_cols].values,\n    full_train_df['opacity']\n)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e5e00dd53d9b6b8c39eb99298dff4bddb0bf24cb"},"cell_type":"markdown","source":"# Load Real Competition Data\nHere we load an image from the real competition and see how our model does"},{"metadata":{"trusted":true,"_uuid":"e04f2cb89ebe60ce5868f681ab86b65e3dbb9100"},"cell_type":"code","source":"comp_dir = '../input/rsna-pneumonia-detection-challenge/'\ncomp_img_dir = os.path.join(comp_dir, 'stage_2_train_images')\nrc_train_df = pd.read_csv(os.path.join(comp_dir, 'stage_2_train_labels.csv'))\nrc_train_df['path'] = rc_train_df['patientId'].map(lambda x: os.path.join(comp_img_dir, '{}.dcm'.format(x)))\nrc_train_df.sample(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8ab5e819b96952d14e03ca951b319c6924be995d"},"cell_type":"code","source":"sort_width_df = rc_train_df.\\\n    groupby('patientId').\\\n    agg({'width': 'sum'}).\\\n    reset_index().\\\n    sort_values('width', ascending=False)\nbig_box = sort_width_df.iloc[0, 0] # first (biggest total width)\nno_box = sort_width_df.iloc[-1, 0] # last (smallest total width)\npneu_pat_df = rc_train_df[rc_train_df['patientId']==big_box]\nhealthy_pat_df = rc_train_df[rc_train_df['patientId']==no_box]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8e92ed1c3fbdabfeab6c6ad50efaef2c7bdfe9dd"},"cell_type":"code","source":"fig, m_axs = plt.subplots(1, 2, figsize = (20, 10))\nsample_df = pd.concat([pneu_pat_df, healthy_pat_df])\nfor c_ax, (c_path, c_rows) in zip(m_axs.flatten(),\n                    sample_df.groupby(['path'])):\n    c_dicom = pydicom.read_file(c_path)\n    c_ax.imshow(c_dicom.pixel_array, cmap='bone')\n    c_ax.set_title('{Target}'.format(**c_rows.iloc[0,:]))\n    for i, (_, c_row) in enumerate(c_rows.dropna().iterrows()):\n        c_ax.plot(c_row['x'], c_row['y'], 's', label='{Target}'.format(**c_row))\n        c_ax.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                                width=c_row['width'],\n                                height=c_row['height'], \n                                 alpha = 0.5))\n        if i==0: c_ax.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"287a20289cb9f0c92f52cab45c054c92b7a596ff"},"cell_type":"code","source":"step_x = full_train_stack.shape[1]*2\nstep_y = full_train_stack.shape[2]*2\noverlap = False\nif overlap:\n    x_steps = np.arange(0, 1024, step_x/2, dtype=int)[1:]\n    y_steps = np.arange(0, 1024, step_y/2, dtype=int)[1:]\nelse:\n    x_steps = np.arange(0, 1024, step_x, dtype=int)+step_x//2\n    y_steps = np.arange(0, 1024, step_y, dtype=int)+step_x//2\nfrom scipy.ndimage import zoom\nopacity_prior_map = zoom(imread(os.path.join(data_dir, 'opacity_prior.tif')), \n                         (training_scale_factor, training_scale_factor),\n                        order=0)\ndef dicom_to_tiles(c_path):\n    c_dicom = pydicom.read_file(c_path)\n    img_rows = []\n    tag_dict = {c_tag: getattr(c_dicom, c_tag, '') for c_tag in DCM_TAG_LIST}\n    for x in x_steps:\n        for y in y_steps:\n            img_rows += [{'x': x, 'y': y, \n                          'img': c_dicom.pixel_array[y-step_y//2:y+step_y//2, x-step_x//2:x+step_x//2],\n                          'opacity_prior': np.mean(opacity_prior_map[y-step_y//2:y+step_y//2, x-step_x//2:x+step_x//2]),\n                         **tag_dict}]\n    out_df = categories_to_indicators(pd.DataFrame(img_rows))\n    # add the intensity channels\n    out_df['Mean_Intensity'] = out_df['img'].map(np.mean)\n    out_df['Std_Intensity'] = out_df['img'].map(np.std)\n    out_df['prediction'] = fitted_pipeline.predict_proba(out_df[model_cols])[:,1]\n    return out_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3e0682620a70be2ffd3376fc0bb3569c61edc7d0"},"cell_type":"code","source":"PREDICTION_CUTOFF = 0.5 # how good a prediction needs to be for us to keep it","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0fa4dad28e480272fcd3a97191562590a864567e"},"cell_type":"code","source":"test_img_path = pneu_pat_df.iloc[0,-1]\ntest_img_tiles_df = dicom_to_tiles(test_img_path)\nprint(test_img_tiles_df.shape[0], 'predictions')\ntest_img_tiles_df.sample(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1065963a5a13a1af6f7ee42490f0fcda2c04645d"},"cell_type":"code","source":"from matplotlib import animation, rc\nrc('animation', html='jshtml', embed_limit=100)\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize = (10, 5))\nc_dicom = pydicom.read_file(test_img_path)\nax1.imshow(c_dicom.pixel_array, cmap='bone')\nax1.set_title('Ground Truth')\nfor i, (_, c_row) in enumerate(pneu_pat_df.dropna().iterrows()):\n    ax1.plot(c_row['x'], c_row['y'], 's', label='{Target}'.format(**c_row))\n    ax1.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                            width=c_row['width'],\n                            height=c_row['height'], \n                             alpha = 0.5))\n    if i==0: ax1.legend()\n\nax2.imshow(c_dicom.pixel_array, cmap='bone')\nax2.set_title('Prediction')\nax2.set_xlim(0, 1024)\nax2.set_ylim(1024, 0)\npoint_list = list(test_img_tiles_df.sort_values(['x', 'y']).dropna().iterrows())\n\ndef add_point(i):\n    _, c_row = point_list[i]\n    c_color = plt.cm.hot(np.array([c_row['prediction']]))\n    artists = [ax2.scatter(c_row['x'], \n                        c_row['y'],\n                           c=c_color\n                       )]\n    artists += [ax2.add_patch(Rectangle(\n        xy=(c_row['x']-step_x//2,\n            c_row['y']-step_y//2),\n        width=step_x,\n        height=step_y, \n        color='red' if c_row['prediction']>PREDICTION_CUTOFF else 'blue',\n        alpha=c_row['prediction'] if c_row['prediction']>PREDICTION_CUTOFF else 0.1))]\n    return tuple(artists)\nani = animation.FuncAnimation(fig, add_point, range(len(point_list)), interval=50)\nani","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ef510639d8db780386a859b00cedc4e587c86615"},"cell_type":"code","source":"from IPython.display import HTML\n# reset the axes\nax2.cla()\nax2.imshow(c_dicom.pixel_array, cmap='bone')\nax2.set_xlim(0, 1024)\nax2.set_ylim(1024, 0)\nax2.set_title('Prediction')\nani.save('step_by_step_prediction.gif', writer='imagemagick')\nHTML('<img src=\"step_by_step_prediction.gif\"/>')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"04fd864db0e55847163ee533f9835b4b7ebe5155"},"cell_type":"markdown","source":"# Try the Healthy Patient"},{"metadata":{"trusted":true,"_uuid":"cd0901a548e0f45715ee94306652189415031a1f"},"cell_type":"code","source":"test_img_path = healthy_pat_df.iloc[0,-1]\ntest_img_tiles_df = dicom_to_tiles(test_img_path)\ntest_img_tiles_df.sample(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d839ea57fb78218231b115a678ed601587acb558"},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (20, 10))\nc_dicom = pydicom.read_file(test_img_path)\nax1.imshow(c_dicom.pixel_array, cmap='bone')\nax1.set_title('Ground Truth')\nfor i, (_, c_row) in enumerate(healthy_pat_df.dropna().iterrows()):\n    ax1.plot(c_row['x'], c_row['y'], 's', label='{Target}'.format(**c_row))\n    ax1.add_patch(Rectangle(xy=(c_row['x'], c_row['y']),\n                            width=c_row['width'],\n                            height=c_row['height'], \n                             alpha = 0.5))\n    if i==0: ax1.legend()\n\nax2.imshow(c_dicom.pixel_array, cmap='bone')\nax2.set_title('Prediction')\nfor i, (_, c_row) in enumerate(test_img_tiles_df.dropna().iterrows()):\n    ax2.plot(c_row['x'], c_row['y'], 's', label='{prediction:.1%}'.format(**c_row))\n    if c_row['prediction']>0.75:\n        ax2.add_patch(Rectangle(xy=(c_row['x']-step_x//2, \n                                    c_row['y']-step_y//2),\n                                width=step_x,\n                                height=step_y, \n                                 alpha = c_row['prediction']/2))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2a00a2bb4d5323e59afec806ab3dc16ef0a9fd22"},"cell_type":"markdown","source":"## Everything seems good\nThe model seems to work, we get many more predictions for the sick patient than the healthy patient we can now try to run all of the test data and see what happens"},{"metadata":{"_uuid":"9f4009cc3949f213beaadbd0c3da424cb951291d"},"cell_type":"markdown","source":"# Running everything\nWe can now run the predictions on all of the images"},{"metadata":{"trusted":true,"_uuid":"60a9aba9bd5e3e13a461ac5da1b10071679280ad"},"cell_type":"code","source":"from glob import glob\nfrom tqdm import tqdm_notebook\nall_test_images = glob(os.path.join(comp_dir, 'stage_2_test_images', '*.dcm'))\nall_predictions_list = {}\nfor c_path in tqdm_notebook(all_test_images):\n    patient_id = os.path.splitext(os.path.basename(c_path))[0]\n    # keep only good predictions and remove the image column\n    all_predictions_list[patient_id] = dicom_to_tiles(c_path).\\\n        query('prediction>{}'.format(PREDICTION_CUTOFF)).\\\n        drop(['img'],1) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"018d0fb13eda662905892de54bfadbca5f696617"},"cell_type":"code","source":"def pkg_result_string(in_df):\n    out_list = in_df[['prediction', 'x', 'y']].apply(lambda c_row: \n                                                     '%2.2f %d %d %d %d' % (c_row['prediction'], \n                                                                          c_row['x']-step_x//2,\n                                                                          c_row['y']-step_y//2,\n                                                                          step_x, \n                                                                          step_y), 1).values.tolist()\n    return ' '.join(out_list)\n#pkg_result_string(all_predictions_list['2d5d8ecc-3ee3-4c9b-bb86-3f614a079585'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"464c3e7bd2446cdb29985f21ae51be3bc08d0a2e"},"cell_type":"code","source":"submission_df = pd.DataFrame({'patientId': list(all_predictions_list.keys())})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"877d824786c39452fd160ab36d87bfdac17c0711"},"cell_type":"code","source":"submission_df['PredictionString'] = submission_df['patientId'].\\\n    map(lambda x: pkg_result_string(all_predictions_list[x]))\nsubmission_df.sample(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a38f3e496d74b55670087c3ee57cd7ec20b5a24"},"cell_type":"code","source":"submission_df.to_csv('submission_rsna_competition.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"048a66fc7ad224cc41a93be5500ed822211cea26"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}