{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nimport pydicom\nimport cv2\nfrom tqdm import tqdm\nfrom tensorflow.keras.models import load_model\nimport glob\n\nfrom tensorflow.keras.applications import InceptionResNetV2\nfrom tensorflow.keras.layers import Dense, Dropout, GlobalAveragePooling2D\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import time \n\nstart_time = time.time()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!cp ../input/gdcm-conda-install/gdcm.tar .\n\n!tar -xvzf gdcm.tar\n!conda install --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2\nprint(\"done\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"IMAGE_SIZE = 512\nMODELS = dict()\nPATH_DATA = '../input/rsna-str-pulmonary-embolism-detection/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(dcm):\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = dcm.pixel_array.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    intercept = dcm.RescaleIntercept\n    slope = dcm.RescaleSlope\n        \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n            \n    image += np.int16(intercept)\n    \n    del intercept, slope, dcm\n    \n    return np.array(image, dtype=np.int16)\n\ndef set_lungwin(img, hu=[-100, 550]):\n    lungwin = np.array(hu)\n    newimg = (img-lungwin[0]) / (lungwin[1]-lungwin[0])\n    newimg[newimg < 0] = 0\n    newimg[newimg > 1] = 1\n    newimg = (newimg * 255).astype('uint8')\n    \n    del img \n    \n    return newimg\n\ndef get_img(d):\n    #d = pydicom.dcmread(path)\n    im = set_lungwin(get_pixels_hu(d))\n    im = cv2.resize(im, (IMAGE_SIZE, IMAGE_SIZE))\n    im = cv2.cvtColor(im,cv2.COLOR_GRAY2RGB)\n    print()\n    del d \n    \n    return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_models():\n    print('load models')\n    MODELS['pe_present'] = load_model(\"../input/models-v1/pe_present_img_inception.h5\")\n    MODELS['ratio'] = load_model(\"../input/models-v1/TEP_ratio_model_v2.h5\")\n    MODELS['side'] = load_model(\"../input/models-v1/TEP_side_model_v2.h5\")\n    MODELS['chronic'] = load_model(\"../input/models-v1/TEP_chronic_model_v1.h5\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_slices(exam_id, series_id):\n    path_patient = PATH_DATA + 'test/' + exam_id + \"/\" + series_id + \"/\"\n    dcm_path = glob.glob(path_patient + \"*\")\n    slices = [pydicom.dcmread(file) for file in dcm_path]\n    \n    del dcm_path, path_patient\n    \n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    \n    start_i = round(len(slices)*0.22)\n    end_i = round(len(slices)*0.85)\n    \n    slices_ini = slices[0:start_i]\n    slices_end = slices[end_i:len(slices)]\n    \n    slices = slices[start_i:end_i]\n    \n    sop_uids = [s.SOPInstanceUID for s in slices]\n    \n    sop_uids_ini = [s.SOPInstanceUID for s in slices_ini]\n    sop_uids_end = [s.SOPInstanceUID for s in slices_end]\n    \n    del slices_ini, slices_end\n    \n    sop_uids = pd.DataFrame(sop_uids, columns=['SOPInstanceUID'])\n    \n    sop_uids['StudyInstanceUID'] = exam_id\n    sop_uids['SeriesInstanceUID'] = series_id\n    \n    sop_uids_ini = pd.DataFrame(sop_uids_ini, columns=['SOPInstanceUID'])\n    sop_uids_ini['StudyInstanceUID'] = exam_id\n    sop_uids_ini['SeriesInstanceUID'] = series_id\n    \n    sop_uids_end = pd.DataFrame(sop_uids_end, columns=['SOPInstanceUID'])\n    sop_uids_end['StudyInstanceUID'] = exam_id\n    sop_uids_end['SeriesInstanceUID'] = series_id\n    \n    del exam_id, series_id\n    \n    return sop_uids, sop_uids_ini, sop_uids_end, slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_predict_columns(df):\n    df[\"pe_present_on_image_pred\"] = 0\n    \n    df[\"rv_lv_ratio_gte_1_pred\"] = 0\n    df[\"rv_lv_ratio_lt_1_pred\"] = 0\n    \n    df[\"leftsided_pe_pred\"] = 0\n    df[\"rightsided_pe_pred\"] = 0\n    df[\"central_pe_pred\"] = 0\n\n    df[\"acute_pe_pred\"] = 0\n    df[\"chronic_pe_pred\"] = 0\n    df[\"acute_and_chronic_pe_pred\"] = 0\n    \n    return df\n\ndef predict_dataset(df, model_pe_present, model_ratio, model_side, model_chronic, slices):\n    \n    df = create_predict_columns(df)\n    \n    N_STEP = len(df)//25\n    B_STEP = len(df)%25\n    START_STEP = B_STEP//2\n    \n    ct = 0\n    pred_pe_present = 0.00001\n    pred_ratio, pred_side, pred_chronic = 0.00001, 0.00001, 0.00001\n    instance, series = 0, 0\n    for idx in df.index:\n        p = PATH_DATA + 'test/' + df.loc[idx,'StudyInstanceUID'] + '/' + df.loc[idx,'SeriesInstanceUID'] + '/' + df.loc[idx,'SOPInstanceUID'] + '.dcm'\n        \n        del instance, series\n        \n        instance = df.loc[idx,'SOPInstanceUID']\n        series = df.loc[idx,'SeriesInstanceUID']\n        \n        if ct==0 or (ct-B_STEP)%N_STEP==0:\n            images = []\n            s = slices[ct]\n            im = get_img(s)\n            images = np.array([im])\n            \n            ct+=1\n           \n            del s, im, pred_pe_present\n            \n            pred_pe_present = model_pe_present.predict(images)\n            \n            \n            df.loc[idx, \"pe_present_on_image_pred\"] = pred_pe_present[0][0]\n            sub.loc[sub[\"id\"]==instance, \"label\"] = pred_pe_present[0][0]\n\n            if pred_pe_present[0][0] > TH_POS_IMG:\n                \n                del pred_ratio, pred_side, pred_chronic\n                \n                pred_ratio  = model_ratio.predict(images)\n                pred_side = model_side.predict(images)\n                pred_chronic = model_chronic.predict(images)\n                \n                mean_negative_exam = (pred_ratio[0][2] + pred_side[0][3] + pred_chronic[0][3])/3\n                \n                if mean_negative_exam < TH_POS_IMG:\n                    pred_pe_present[0][0] = min(mean_negative_exam, 0.49)\n                    sub.loc[sub[\"id\"]==instance, \"label\"] = pred_pe_present[0][0]\n                    continue\n                \n                df.loc[idx, \"rv_lv_ratio_gte_1_pred\"] = pred_ratio[0][0] \n                df.loc[idx, \"rv_lv_ratio_lt_1_pred\"] = pred_ratio[0][1]\n\n                df.loc[idx, \"leftsided_pe_pred\"] = pred_side[0][0]\n                df.loc[idx, \"rightsided_pe_pred\"] = pred_side[0][1]\n                df.loc[idx, \"central_pe_pred\"] = pred_side[0][2]\n\n                df.loc[idx, \"acute_pe_pred\"] = pred_chronic[0][0]\n                df.loc[idx, \"chronic_pe_pred\"] = pred_chronic[0][1]\n                df.loc[idx, \"acute_and_chronic_pe_pred\"] = pred_chronic[0][2]\n                \n                del images\n            \n            else:\n                del images\n        \n\n        else:\n            ct+=1\n            df.loc[idx, \"pe_present_on_image_pred\"] = pred_pe_present[0][0]\n            sub.loc[sub[\"id\"]==instance, \"label\"] = pred_pe_present[0][0]\n\n            if pred_pe_present[0][0] > TH_POS_IMG:\n                df.loc[idx, \"rv_lv_ratio_gte_1_pred\"] = pred_ratio[0][0]\n                df.loc[idx, \"rv_lv_ratio_lt_1_pred\"] = pred_ratio[0][1]\n\n                df.loc[idx, \"leftsided_pe_pred\"] = pred_side[0][0]\n                df.loc[idx, \"rightsided_pe_pred\"] = pred_side[0][1]\n                df.loc[idx, \"central_pe_pred\"] = pred_side[0][2]\n\n                df.loc[idx, \"acute_pe_pred\"] = pred_chronic[0][0]\n                df.loc[idx, \"chronic_pe_pred\"] = pred_chronic[0][1]\n                df.loc[idx, \"acute_and_chronic_pe_pred\"] = pred_chronic[0][2]\n    \n    del slices\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"load_models()\nprint('load data')\n\nfrom os import path\nif path.exists('../input/rsna-str-pulmonary-embolism-detection/train'):\n    test=pd.read_csv(PATH_DATA+'/test.csv').head(2000)#\nelse:\n    test=pd.read_csv(PATH_DATA+'/test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#sub.label = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def treat_pred(exam_pred, prob_neg):\n    if prob_neg < 0.5: #exam positive\n        \n        #Let's solve rule1a\n        if exam_pred['rv_lv_ratio_lt_1_pred']  >  0.5  and exam_pred['rv_lv_ratio_gte_1_pred'] >  0.5:\n            # The greater keeps the same, the lower will be the opposite of the greater\n            # This case should not accure in our case \n            if exam_pred['rv_lv_ratio_lt_1_pred'] > exam_pred['rv_lv_ratio_gte_1_pred']:\n                exam_pred['rv_lv_ratio_gte_1_pred'] = 1 - exam_pred['rv_lv_ratio_lt_1_pred']\n            else:\n                exam_pred['rv_lv_ratio_lt_1_pred'] = 1 - exam_pred['rv_lv_ratio_gte_1_pred']\n                \n        elif exam_pred['rv_lv_ratio_lt_1_pred']  <= 0.5  and exam_pred['rv_lv_ratio_gte_1_pred'] <= 0.5:\n            # This case might happen too, it is when the model of ratio thinks that the exam is negative\n            # To solve this we will take the biggest pred and make the max with prob_neg\n            # This solution is lame but the model too, so we cannot relie on it, yet \n            if exam_pred['rv_lv_ratio_lt_1_pred'] > exam_pred['rv_lv_ratio_gte_1_pred']:\n                exam_pred['rv_lv_ratio_lt_1_pred'] = 1 - prob_neg\n                exam_pred['rv_lv_ratio_gte_1_pred'] = prob_neg\n            else:\n                exam_pred['rv_lv_ratio_gte_1_pred'] = 1 - prob_neg\n                exam_pred['rv_lv_ratio_lt_1_pred'] = prob_neg\n        \n        #Let's solve rule1b\n        if exam_pred['central_pe_pred'] <= 0.5 and exam_pred['rightsided_pe_pred'] <= 0.5 and exam_pred['leftsided_pe_pred'] <= 0.5:\n            # We use the same idea, maybe we can improve this because the model is really better\n            max_pred = max(exam_pred['central_pe_pred'], exam_pred['rightsided_pe_pred'], exam_pred['leftsided_pe_pred'])\n            for c in side:\n                if exam_pred[c] == max_pred:\n                    exam_pred[c] = 1 - prob_neg\n        \n        # Let's solve rule1c\n        if exam_pred['acute_and_chronic_pe_pred'] > 0.5 and exam_pred['chronic_pe_pred'] > 0.5:\n            # Really poor chances that this case happens since our model has got a huge proportion to predict acute_pe, but we if it happens we just gonna keep the max\n            # since if the model gives a high pred for these two classes it would mean that it is REALLY NOT acute\n            if exam_pred['acute_and_chronic_pe_pred'] > exam_pred['chronic_pe_pred']:\n                exam_pred['chronic_pe_pred'] = 1 - exam_pred['acute_and_chronic_pe_pred']\n            else:\n                exam_pred['acute_and_chronic_pe_pred'] = 1 - exam_pred['chronic_pe_pred']\n                \n        # Since we give 0 everytime for indeterminate we have no risk to break rule 1d\n        \n    else: # Exam negative\n        # It seems to me that the rule 2a cannot be broken but I will keep an eye on it \n        # Rule 2b will be solved by putting a filter, every pred superior to the probability of being positive will be equal to the prob of being positive\n        # This rule makes total sense and can be kept for long term\n        for c in classes:\n            exam_pred[c] = min(exam_pred[c], 1 - prob_neg)\n    \n    return exam_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import gc \n\nTH_POS_IMG = 0.5\nTH_POS = 0.90\n\n# Para cada paciente:\nfor p in tqdm(test['StudyInstanceUID'].unique()):\n    \n    if (time.time() - start_time) > 8.5*60*60:\n        break\n\n    gc.collect()\n    \n    # pega slices ordenados, cortando % inicial e final\n    sop_uids, sop_uids_ini, sop_uids_end, slices = get_slices(p, test[test['StudyInstanceUID']==p]['SeriesInstanceUID'].values[0])\n\n    for idx in sop_uids_ini.SOPInstanceUID:\n        sub.loc[sub[\"id\"]==idx, \"label\"] = 0.000001\n    \n    for idx in sop_uids_end.SOPInstanceUID:\n        sub.loc[sub[\"id\"]==idx, \"label\"] = 0.000001\n    \n    del sop_uids_ini, sop_uids_end\n    \n    # predição dos slices desse exame:\n    exam_predictions = predict_dataset(sop_uids, MODELS['pe_present'], MODELS['ratio'], MODELS['side'], MODELS['chronic'], slices)\n    \n    del sop_uids\n    \n    classes = ['rv_lv_ratio_gte_1_pred', 'rv_lv_ratio_lt_1_pred', 'leftsided_pe_pred', 'rightsided_pe_pred', 'central_pe_pred', \n               'chronic_pe_pred', 'acute_and_chronic_pe_pred']\n    \n    ratio = ['rv_lv_ratio_gte_1_pred', 'rv_lv_ratio_lt_1_pred']\n    \n    side = ['leftsided_pe_pred', 'rightsided_pe_pred', 'central_pe_pred']\n    \n    chronic = ['chronic_pe_pred', 'acute_and_chronic_pe_pred']\n    \n    max_score = exam_predictions.pe_present_on_image_pred.max()\n\n    prob_neg = 1 - max_score\n    \n    del max_score \n    \n    idx = p + '_negative_exam_for_pe'\n    sub.loc[sub[\"id\"]==idx, \"label\"] = prob_neg\n    \n    exam_pred = {}\n    for c in classes:\n        pr = exam_predictions[exam_predictions[c]>0][c].mean()\n        exam_pred[c] = pr\n    \n    exam_pred_treated = treat_pred(exam_pred, prob_neg)\n    \n    del exam_pred, prob_neg\n    \n    for c in classes:\n        idx = p + '_' + c.replace('_pred','')\n        sub.loc[sub[\"id\"]==idx, \"label\"] = exam_pred_treated[c]\n    \n    del exam_pred_treated\n    \n    idx = p + '_' + 'indeterminate'.replace('_pred','')\n    \n    sub.loc[sub[\"id\"]==idx, \"label\"] = 0.020484822355039723\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def check_consistency(sub, test):\n    \n    '''\n    Checks label consistency and returns the errors\n    \n    Args:\n    sub   = submission dataframe (pandas)\n    test  = test.csv dataframe (pandas)\n    '''\n    \n    # EXAM LEVEL\n    for i in test['StudyInstanceUID'].unique():\n        df_tmp = sub.loc[sub.id.str.contains(i, regex = False)].reset_index(drop = True)\n        df_tmp['StudyInstanceUID'] = df_tmp['id'].str.split('_').str[0]\n        df_tmp['label_type']       = df_tmp['id'].str.split('_').str[1:].apply(lambda x: '_'.join(x))\n        del df_tmp['id']\n        if i == test['StudyInstanceUID'].unique()[0]:\n            df = df_tmp.copy()\n        else:\n            df = pd.concat([df, df_tmp], axis = 0)\n    df_exam = df.pivot(index = 'StudyInstanceUID', columns = 'label_type', values = 'label')\n    \n    # IMAGE LEVEL\n    df_image = sub.loc[sub.id.isin(test.SOPInstanceUID)].reset_index(drop = True)\n    df_image = df_image.merge(test, how = 'left', left_on = 'id', right_on = 'SOPInstanceUID')\n    df_image.rename(columns = {\"label\": \"pe_present_on_image\"}, inplace = True)\n    del df_image['id']\n    \n    # MERGER\n    df = df_exam.merge(df_image, how = 'left', on = 'StudyInstanceUID')\n    ids    = ['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID']\n    labels = [c for c in df.columns if c not in ids]\n    df = df[ids + labels]\n    \n    # SPLIT NEGATIVE AND POSITIVE EXAMS\n    df['positive_images_in_exam'] = df['StudyInstanceUID'].map(df.groupby(['StudyInstanceUID']).pe_present_on_image.max())\n    df_pos = df.loc[df.positive_images_in_exam >  0.5]\n    df_neg = df.loc[df.positive_images_in_exam <= 0.5]\n    \n    # CHECKING CONSISTENCY OF POSITIVE EXAM LABELS\n    rule1a = df_pos.loc[((df_pos.rv_lv_ratio_lt_1  >  0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 >  0.5)) | \n                        ((df_pos.rv_lv_ratio_lt_1  <= 0.5)  & \n                         (df_pos.rv_lv_ratio_gte_1 <= 0.5))].reset_index(drop = True)\n    rule1a['broken_rule'] = '1a'\n    rule1b = df_pos.loc[(df_pos.central_pe    <= 0.5) & \n                        (df_pos.rightsided_pe <= 0.5) & \n                        (df_pos.leftsided_pe  <= 0.5)].reset_index(drop = True)\n    rule1b['broken_rule'] = '1b'\n    rule1c = df_pos.loc[(df_pos.acute_and_chronic_pe > 0.5) & \n                        (df_pos.chronic_pe           > 0.5)].reset_index(drop = True)\n    rule1c['broken_rule'] = '1c'\n    rule1d = df_pos.loc[(df_pos.indeterminate        > 0.5) | \n                        (df_pos.negative_exam_for_pe > 0.5)].reset_index(drop = True)\n    rule1d['broken_rule'] = '1d'\n\n    # CHECKING CONSISTENCY OF NEGATIVE EXAM LABELS\n    rule2a = df_neg.loc[((df_neg.indeterminate        >  0.5)  & \n                         (df_neg.negative_exam_for_pe >  0.5)) | \n                        ((df_neg.indeterminate        <= 0.5)  & \n                         (df_neg.negative_exam_for_pe <= 0.5))].reset_index(drop = True)\n    rule2a['broken_rule'] = '2a'\n    rule2b = df_neg.loc[(df_neg.rv_lv_ratio_lt_1     > 0.5) | \n                        (df_neg.rv_lv_ratio_gte_1    > 0.5) |\n                        (df_neg.central_pe           > 0.5) | \n                        (df_neg.rightsided_pe        > 0.5) | \n                        (df_neg.leftsided_pe         > 0.5) |\n                        (df_neg.acute_and_chronic_pe > 0.5) | \n                        (df_neg.chronic_pe           > 0.5)].reset_index(drop = True)\n    rule2b['broken_rule'] = '2b'\n    \n    # MERGING INCONSISTENT PREDICTIONS\n    errors = pd.concat([rule1a, rule1b, rule1c, rule1d, rule2a, rule2b], axis = 0)\n    \n    # OUTPUT\n    print('Found', len(errors), 'inconsistent predictions')\n    return errors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# max_score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# len(check_consistency(sub, test))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# error = check_consistency(sub, test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# error.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(time.time() - start_time)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}