{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nimport pydicom\nfrom scipy import ndimage\nfrom tqdm import tqdm\nimport sys\nimport matplotlib.pyplot as plt\n\nbasepath = Path('../input/rsna-str-pulmonary-embolism-detection/')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath / 'train.csv')\ntrain['dcmpath'] = str(basepath) + '/train' + '/' + train.StudyInstanceUID + '/' + train.SeriesInstanceUID","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"studies = {}\nlabel_names = ['negative_exam_for_pe', 'qa_motion', 'qa_contrast', 'flow_artifact', 'rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1', 'leftsided_pe', 'chronic_pe', 'true_filling_defect_not_pe', 'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate']\n\nN = 2000000\n\nfor i, row in tqdm(enumerate(train.iterrows())):\n    \n    if i >= N: break\n    \n    study_uid = row[1]['StudyInstanceUID']\n    serie_uid = row[1]['SeriesInstanceUID']\n    slice_uid = row[1]['SOPInstanceUID']\n    \n    if study_uid not in studies.keys():\n        studies[study_uid] = {**{label: row[1][label] for label in label_names}, 'series': {}}\n        \n    series = studies[study_uid]['series']\n    \n    if serie_uid not in series.keys():\n        series[serie_uid] = {'slices': {}}\n        \n    slices = series[serie_uid]['slices']\n    slices[slice_uid] = {'pe_present_on_image': row[1]['pe_present_on_image']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(studies))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"lengths = []\nfor suid in serie_uids:\n    files = np.array(list(Path(suid).iterdir()))\n    lengths.append(len(files))\nprint(np.mean(lengths))\nprint(np.std(lengths))\nprint(np.min(lengths))\nprint(np.max(lengths))\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_scans(study_uids, size=128, depth=64):\n    scans = []\n    #labels = []\n    targets = []\n    for study_uid in study_uids:\n        #labels.append({label: studies[study_uid][label] for label in label_names})\n        serie_uid = list(studies[study_uid]['series'].keys())[0]\n        serie_path = '../input/rsna-str-pulmonary-embolism-detection/train/{0}/{1}'.format(study_uid, serie_uid)\n        slice_paths = list(Path(serie_path).iterdir())\n        #labels[-1]['num_slices'] = len(slice_paths)\n        slice_dcms = np.array([pydicom.dcmread(str(path)) for path in slice_paths])\n        pe_present_on_image = np.array([studies[study_uid]['series'][serie_uid]['slices'][os.path.basename(path).split('.')[0]]['pe_present_on_image'] for path in slice_paths])\n        sortkey = np.argsort([float(s.ImagePositionPatient[2]) for s in slice_dcms])\n        slice_dcms = slice_dcms[sortkey]\n        pe_present_on_image = pe_present_on_image[sortkey]\n        pe_present_on_image_resized = ndimage.interpolation.zoom(pe_present_on_image, 2772./len(pe_present_on_image), order=2)\n        #labels[-1]['pe_present_on_image'] = pe_present_on_image\n        #labels[-1]['pe_present_on_image_resized'] = pe_present_on_image_resized\n        scan = []\n        for dcm in slice_dcms:\n            try:\n                scan.append(dcm.pixel_array.astype(np.float32)*dcm.RescaleSlope + dcm.RescaleIntercept)\n            except Exception as e:\n                print(repr(e))\n            del dcm\n        del slice_dcms\n        scan = np.stack(scan, axis=0)\n        sf = size / 512.\n        df = depth / float(scan.shape[0])\n        scan = ndimage.interpolation.zoom(scan, (df, sf, sf), order=1)\n        scan = np.clip(scan, -1000, 1000)/1000\n        scans.append(scan)\n        del scan\n        targets.append([studies[study_uid]['negative_exam_for_pe'],\n                       studies[study_uid]['indeterminate'],\n                       studies[study_uid]['chronic_pe'],\n                       studies[study_uid]['acute_and_chronic_pe'],\n                       studies[study_uid]['central_pe'],\n                       studies[study_uid]['leftsided_pe'],\n                       studies[study_uid]['rightsided_pe'],\n                       studies[study_uid]['rv_lv_ratio_gte_1'],\n                       studies[study_uid]['rv_lv_ratio_lt_1'],\n                        *pe_present_on_image_resized])\n    targets = np.array(targets)\n    scans = np.array(scans)[..., None]\n    return scans, targets\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls -l\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n = 1200\nwhile n <= 7279:\n    m = n + 100\n    m = min(m, 7279)\n    print('processing {0} to {1}.npy'.format(n, m))\n    scans, targets = get_scans(list(studies.keys())[n:m])\n    np.save('scans_{0}_to_{1}.npy'.format(n, m), scans)\n    np.save('targets_{0}_to_{1}.npy'.format(n, m), targets)\n    del scans\n    del targets\n    if n == 7279: break\n    n = m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N_SCANS = 100\nPERCENT_TRAIN = 70\n\nscans, targets = get_scans(list(studies.keys())[0:0+N_SCANS])\n\nn = int(np.round(PERCENT_TRAIN/100. * scans.shape[0]))\nx_train = scans[:n]\ny_train = targets[:n]\nx_val = scans[n:]\ny_val = targets[n:]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":false,"trusted":true,"collapsed":true},"cell_type":"code","source":"scan = scans[0]\n\nplt.hist(scan.flat, bins=100)\nplt.show()\n\nfig, ax = plt.subplots(5, 4, figsize=(20,20))\nax = ax.flatten()\nfor m in range(20):\n    ax[m].imshow(scan[m], cmap='Blues_r')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow.keras.layers import BatchNormalization, Conv3D, MaxPooling3D, AveragePooling3D, Lambda, UpSampling1D, Conv1D, GlobalMaxPooling3D, Dense ,Flatten, Activation","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SIZE = 128\nDEPTH = 64\n\ntf.keras.backend.clear_session()\n\ninput_tensor = tf.keras.layers.Input(shape=(DEPTH, SIZE, SIZE, 1))\ntarget_input = tf.keras.layers.Input(shape=(2781,))\n\nx = Conv3D(filters=16, kernel_size=(3, 3, 3), padding='valid', strides=(2, 2, 2))(input_tensor)\nx = BatchNormalization()(x)\nx = Activation('relu')(x)\n\nx = Conv3D(filters=24, kernel_size=(3, 3, 3), padding='valid', strides=(1, 2, 2))(x)\nx = BatchNormalization()(x)\nx = Activation('relu')(x)\n\nx = Conv3D(filters=40, kernel_size=(3, 3, 3), padding='valid', strides=(1, 2, 2))(x)\nx = BatchNormalization()(x)\nx = Activation('relu')(x)\n\nx = Conv3D(filters=64, kernel_size=(3, 3, 3), padding='valid', strides=(1, 2, 2))(x)\nx = BatchNormalization()(x)\nx = Activation('relu')(x)\n\nx = Conv3D(filters=128, kernel_size=(3, 3, 3), padding='valid', strides=(1, 2, 2))(x)\nx = BatchNormalization()(x)\nx = Activation('relu')(x)\n\ny = Lambda(lambda t: tf.keras.backend.mean(t, axis=[2, 3]))(x)\ny = UpSampling1D(size=4)(y)\ny = Conv1D(filters=64, kernel_size=5, activation='relu')(y)\ny = UpSampling1D(size=4)(y)\ny = Conv1D(filters=32, kernel_size=5, activation='relu')(y)\ny = UpSampling1D(size=4)(y)\ny = Conv1D(filters=16, kernel_size=5, activation='relu')(y)\ny = UpSampling1D(size=2)(y)\ny = Conv1D(filters=1, kernel_size=5)(y)\ny = Flatten()(y)\n\nz = Conv3D(filters=128, kernel_size=(3, 1, 1), padding='valid', strides=(2, 1, 1))(x)\nz = BatchNormalization()(z)\nz = Activation('relu')(z)\n\nz = GlobalMaxPooling3D()(z)\nz = Dense(9)(z)\n\nclass Loss(tf.keras.layers.Layer):\n    def __init__(self, *args, **kwargs):\n        super(Loss, self).__init__(*args, **kwargs)\n    def call(self, x):\n        \"\"\"\n        x: of shape (batch_size, 9 + size of image vector)\n        \"\"\"\n        \n        gt, pred = x\n        \n        exam_gt = gt[..., :9]\n        exam_pred = pred[..., :9]\n        im_gt = gt[..., 9:]\n        im_pred = pred[..., 9:]\n        \n        weights = tf.constant([0.0736196319, 0.09202453988, 0.1042944785, 0.1042944785, 0.1877300613, 0.06257668712, 0.06257668712, 0.2346625767, 0.0782208589])\n        weights = tf.constant([1, 0, 0, 0, 0, 0, 0, 0, 0], tf.float32)\n        weights = weights[None, ...]\n        \n        exam_loss = tf.reduce_sum(weights * tf.nn.sigmoid_cross_entropy_with_logits(labels=exam_gt, logits=exam_pred), axis=-1)\n        \n        #tf.print(exam_gt[0])\n        #tf.print(tf.nn.sigmoid(exam_pred[0]))\n        #tf.print(exam_loss[0])\n        \n        #q = tf.reduce_mean(im_gt, axis=-1)\n        #image_loss = 0.25 * 0.07361963 * q * tf.reduce_sum(tf.nn.sigmoid_cross_entropy_with_logits(labels=im_gt, logits=im_pred), axis=-1)\n        image_loss = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(labels=im_gt, logits=im_pred), axis=-1)\n        \n        loss = exam_loss + image_loss\n        \n        #loss = loss / tf.reduce_sum(weights) + 0.07361963 * q\n        loss = tf.reduce_mean(loss) # mean across batch\n        \n        self.add_loss(loss)\n        self.add_metric(exam_loss, aggregation='mean', name='exam_loss')\n        self.add_metric(image_loss, aggregation='mean', name='image_loss')\n        \n        return loss\n\npred = Lambda(lambda x: tf.keras.backend.concatenate([x[0], x[1]], axis=-1))([z, y])\nloss = Loss()([target_input, pred])\n\nmodel = tf.keras.Model(inputs=[input_tensor, target_input], outputs=[pred, loss])\n#model.summary()\nmodel.compile(optimizer=tf.keras.optimizers.Adam(1e-3))\n\nmodel.fit(x=[x_train, y_train], y=None, batch_size=16, epochs=20, validation_data=[x_val, y_val], validation_freq=2)\n\n#print(model(scans))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.compile(optimizer=tf.keras.optimizers.Adam(1e-4))\nmodel.fit(x=[x_train, y_train], y=None, batch_size=16, epochs=20, validation_data=[x_val, y_val], validation_freq=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.save_weights('model.h5')\nx = 5\npred, loss = model.predict([scans[x:x+1], targets[x:x+1]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nx = 2\npred, loss = model.predict([scans[x:x+1], targets[x:x+1]])\n\npred_sig = 1/(1 + np.exp(-pred)) \n\nplt.plot(pred_sig[0])\nplt.plot(targets[x])\nplt.show()\n\nplt.plot(pred_sig[0, :9])\nplt.plot(targets[x, :9])\nplt.show()\n\nwith np.printoptions(precision=3, suppress=True):\n    print(np.concatenate([pred_sig[..., :9], targets[x:x+1][..., :9]], axis=0))\n","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}