{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np # linear algebra\nimport tensorflow as tf # for tensorflow based registration\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom skimage.util.montage import montage2d\nimport os\nfrom cv2 import imread, createCLAHE # read and equalize images\nimport cv2\nfrom glob import glob\n%matplotlib inline\nimport matplotlib.pyplot as plt","metadata":{"_cell_guid":"f42f6560-edf0-4efb-85a6-6e945e50895b","_uuid":"3300a1edbf2e8122d88093998eb503a6fab8a719","execution":{"iopub.status.busy":"2023-04-07T09:43:01.411427Z","iopub.execute_input":"2023-04-07T09:43:01.411936Z","iopub.status.idle":"2023-04-07T09:43:04.068215Z","shell.execute_reply.started":"2023-04-07T09:43:01.411839Z","shell.execute_reply":"2023-04-07T09:43:04.067546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cxr_paths = glob(os.path.join('..', 'input', 'pulmonary-chest-xray-abnormalities',\n                              'Montgomery', 'MontgomerySet', '*', '*.png'))\ncxr_images = [(c_path, \n               [os.path.join('/'.join(c_path.split('/')[:-2]),'ManualMask','leftMask', os.path.basename(c_path)),\n               os.path.join('/'.join(c_path.split('/')[:-2]),'ManualMask','rightMask', os.path.basename(c_path))]\n              ) for c_path in cxr_paths]\nprint('CXR Images', len(cxr_paths), cxr_paths[0])\nprint(cxr_images[0])","metadata":{"_cell_guid":"9a342cdc-0823-490d-9a3a-a53fb7c33727","_uuid":"fe804e7c294e2d290e27b037bf1ba56177abab70","execution":{"iopub.status.busy":"2023-04-07T09:43:04.069240Z","iopub.execute_input":"2023-04-07T09:43:04.069546Z","iopub.status.idle":"2023-04-07T09:43:04.410910Z","shell.execute_reply.started":"2023-04-07T09:43:04.069499Z","shell.execute_reply":"2023-04-07T09:43:04.410135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skimage.io import imread as imread_raw\nfrom skimage.transform import resize\nimport warnings\nfrom tqdm import tqdm\nwarnings.filterwarnings('ignore', category=UserWarning, module='skimage') # skimage is really annoying\nOUT_DIM = (512, 512)\ndef imread(in_path, apply_clahe = False):\n    img_data = imread_raw(in_path)\n    n_img = (255*resize(img_data, OUT_DIM, mode = 'constant')).clip(0,255).astype(np.uint8)\n    if apply_clahe:\n        clahe_tool = createCLAHE(clipLimit=2.0, tileGridSize=(16,16))\n        n_img = clahe_tool.apply(n_img)\n    return np.expand_dims(n_img, -1)","metadata":{"_cell_guid":"8e0abf66-0f34-4289-8eb6-df6ddd3c2f11","_uuid":"90d2424684780ec518c8ef37f7f988cc24f4ad08","execution":{"iopub.status.busy":"2023-04-07T09:43:04.412407Z","iopub.execute_input":"2023-04-07T09:43:04.412795Z","iopub.status.idle":"2023-04-07T09:43:04.693658Z","shell.execute_reply.started":"2023-04-07T09:43:04.412715Z","shell.execute_reply":"2023-04-07T09:43:04.692947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_vol, seg_vol = [], []\nfor img_path, s_paths in tqdm(cxr_images):\n    img_vol += [imread(img_path)]    \n    seg_vol += [np.max(np.stack([imread(s_path, apply_clahe = False) for s_path in s_paths],0),0)]\nimg_vol = np.stack(img_vol,0)\nseg_vol = np.stack(seg_vol,0)\nprint('Images', img_vol.shape, 'Segmentations', seg_vol.shape)","metadata":{"_uuid":"284ee84254f9ee082341d4f6125be053fa89dd0e","execution":{"iopub.status.busy":"2023-04-07T09:43:04.694878Z","iopub.execute_input":"2023-04-07T09:43:04.695358Z","iopub.status.idle":"2023-04-07T09:46:17.749847Z","shell.execute_reply.started":"2023-04-07T09:43:04.695292Z","shell.execute_reply":"2023-04-07T09:46:17.749094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed(2018)\nt_img, m_img = img_vol[0], seg_vol[0]\n\nfig, (ax_img, ax_mask) = plt.subplots(1,2, figsize = (12, 6))\nax_img.imshow(np.clip(255*t_img, 0, 255).astype(np.uint8) if t_img.shape[2]==3 else t_img[:,:,0],\n              interpolation = 'none', cmap = 'bone')\nax_mask.imshow(m_img[:,:,0], cmap = 'bone')","metadata":{"_cell_guid":"a1367d63-4c7b-4e47-b19f-9a26d1f89476","scrolled":true,"_uuid":"6ed6489bbd618fb419ceca2bbd8c300694f1d4ed","execution":{"iopub.status.busy":"2023-04-07T09:46:17.751275Z","iopub.execute_input":"2023-04-07T09:46:17.752002Z","iopub.status.idle":"2023-04-07T09:46:18.018946Z","shell.execute_reply.started":"2023-04-07T09:46:17.751947Z","shell.execute_reply":"2023-04-07T09:46:18.018206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Make a Simple Model\nHere we make a simple U-Net to create the lung segmentations","metadata":{"_uuid":"40e1de3dca248a362fec0563d582a0536efd44bd"}},{"cell_type":"code","source":"from keras.layers import Conv2D, Activation, Input, UpSampling2D, concatenate, BatchNormalization\nfrom keras.layers import LeakyReLU\nfrom keras.initializers import RandomNormal\ndef c2(x_in, nf, strides=1):\n    x_out = Conv2D(nf, kernel_size=3, padding='same',\n                   kernel_initializer='he_normal', strides=strides)(x_in)\n    x_out = LeakyReLU(0.2)(x_out)\n    return x_out\ndef unet_enc(vol_size, enc_nf, pre_filter = 8):\n    src = Input(shape=vol_size + (1,), name = 'EncoderInput')\n    # down-sample path.\n    x_in = BatchNormalization(name = 'NormalizeInput')(src)\n    x_in = c2(x_in, pre_filter, 1)\n    x0 = c2(x_in, enc_nf[0], 2)  \n    x1 = c2(x0, enc_nf[1], 2)  \n    x2 = c2(x1, enc_nf[2], 2)  \n    x3 = c2(x2, enc_nf[3], 2) \n    return Model(inputs = [src], \n                outputs = [x_in, x0, x1, x2, x3],\n                name = 'UnetEncoder')","metadata":{"_uuid":"3ccfe7e7b7f40caa7fa6aa9ba8d9803538c2bca3","execution":{"iopub.status.busy":"2023-04-07T09:46:18.020240Z","iopub.execute_input":"2023-04-07T09:46:18.020666Z","iopub.status.idle":"2023-04-07T09:46:18.379135Z","shell.execute_reply.started":"2023-04-07T09:46:18.020512Z","shell.execute_reply":"2023-04-07T09:46:18.378386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras.models import Model\nfrom keras import layers\ndef unet(vol_size, enc_nf, dec_nf, full_size=True, edge_crop=48):\n    \"\"\"\n    unet network for voxelmorph \n    Args:\n        vol_size: volume size. e.g. (256, 256, 256)\n        enc_nf: encoder filters. right now it needs to be to 1x4.\n            e.g. [16,32,32,32]\n            TODO: make this flexible.\n        dec_nf: encoder filters. right now it's forced to be 1x7.\n            e.g. [32,32,32,32,8,8,3]\n            TODO: make this flexible.\n        full_size\n    \"\"\"\n\n    # inputs\n    raw_src = Input(shape=vol_size + (1,), name = 'ImageInput')\n    src = layers.GaussianNoise(0.25)(raw_src)\n    enc_model = unet_enc(vol_size, enc_nf)\n    # run the same encoder on the source and the target and concatenate the output at each level\n    x_in, x0, x1, x2, x3 = [s_enc for s_enc in enc_model(src)]\n\n    x = c2(x3, dec_nf[0])\n    x = UpSampling2D()(x)\n    x = concatenate([x, x2])\n    x = c2(x, dec_nf[1])\n    x = UpSampling2D()(x)\n    x = concatenate([x, x1])\n    x = c2(x, dec_nf[2])\n    x = UpSampling2D()(x)\n    x = concatenate([x, x0])\n    x = c2(x, dec_nf[3])\n    x = c2(x, dec_nf[4])\n    x = UpSampling2D()(x)\n    x = concatenate([x, x_in])\n    x = c2(x, dec_nf[5])\n\n    # transform the results into a flow.\n    y_seg = Conv2D(1, kernel_size=3, padding='same', name='lungs', activation='sigmoid')(x)\n    y_seg = layers.Cropping2D((edge_crop, edge_crop))(y_seg)\n    y_seg = layers.ZeroPadding2D((edge_crop, edge_crop))(y_seg)\n    # prepare model\n    model = Model(inputs=[raw_src], outputs=[y_seg])\n    return model","metadata":{"_uuid":"180055fc4b225ea5a39db3bc01085f920b70e139","execution":{"iopub.status.busy":"2023-04-07T09:46:18.380412Z","iopub.execute_input":"2023-04-07T09:46:18.380845Z","iopub.status.idle":"2023-04-07T09:46:18.457673Z","shell.execute_reply.started":"2023-04-07T09:46:18.380796Z","shell.execute_reply":"2023-04-07T09:46:18.456826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# use the predefined depths\nnf_enc=[16,32,32,32]\nnf_dec=[32,32,32,32,32,16,16,2]\nnet = unet(OUT_DIM, nf_enc, nf_dec)\n# ensure the model roughly works\na= net.predict([np.zeros((1,)+OUT_DIM+(1,))])\nprint(a.shape)\nnet.summary()","metadata":{"_uuid":"75391ccba53a3e9beeebd3c74d3a02fc2c5fdd93","execution":{"iopub.status.busy":"2023-04-07T09:46:18.459058Z","iopub.execute_input":"2023-04-07T09:46:18.459555Z","iopub.status.idle":"2023-04-07T09:46:26.154976Z","shell.execute_reply.started":"2023-04-07T09:46:18.459504Z","shell.execute_reply":"2023-04-07T09:46:26.154244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras.optimizers import Adam\nimport keras.backend as K\nfrom keras.optimizers import Adam\nfrom keras.losses import binary_crossentropy\n\nreg_param = 1.0\nlr = 2e-4\ndice_bce_param = 0.0\nuse_dice = True\n\ndef dice_coef(y_true, y_pred, smooth=1):\n    intersection = K.sum(y_true * y_pred, axis=[1,2,3])\n    union = K.sum(y_true, axis=[1,2,3]) + K.sum(y_pred, axis=[1,2,3])\n    return K.mean( (2. * intersection + smooth) / (union + smooth), axis=0)\ndef dice_p_bce(in_gt, in_pred):\n    return dice_bce_param*binary_crossentropy(in_gt, in_pred) - dice_coef(in_gt, in_pred)\ndef true_positive_rate(y_true, y_pred):\n    return K.sum(K.flatten(y_true)*K.flatten(K.round(y_pred)))/K.sum(y_true)\n\nnet.compile(optimizer=Adam(lr=lr), \n              loss=[dice_p_bce], \n           metrics = [true_positive_rate, 'binary_accuracy'])","metadata":{"_uuid":"66a299e9dc4d8618f2a752fe9e9bd7849959c76b","execution":{"iopub.status.busy":"2023-04-07T09:46:26.157218Z","iopub.execute_input":"2023-04-07T09:46:26.157516Z","iopub.status.idle":"2023-04-07T09:46:26.227966Z","shell.execute_reply.started":"2023-04-07T09:46:26.157476Z","shell.execute_reply":"2023-04-07T09:46:26.227362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Training Data Generator\nHere we make a tool to generate training data from the X-ray scans","metadata":{"_uuid":"9d48d2a5f7610c97c374f517cfb23f77aaedc138"}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\ntrain_vol, test_vol, train_seg, test_seg = train_test_split((img_vol-127.0)/127.0, \n                                                            (seg_vol>127).astype(np.float32), \n                                                            test_size = 0.2, \n                                                            random_state = 2018)\nprint('Train', train_vol.shape, 'Test', test_vol.shape, test_vol.mean(), test_vol.max())\nprint('Seg', train_seg.shape, train_seg.max(), np.unique(train_seg.ravel()))\nfig, (ax1, ax1hist, ax2, ax2hist) = plt.subplots(1, 4, figsize = (20, 4))\nax1.imshow(test_vol[0, :, :, 0])\nax1hist.hist(test_vol.ravel())\nax2.imshow(test_seg[0, :, :, 0]>0.5)\nax2hist.hist(train_seg.ravel());","metadata":{"_uuid":"063195f94409876ac24e6cb9f0b667cd9974d289","execution":{"iopub.status.busy":"2023-04-07T09:46:26.229134Z","iopub.execute_input":"2023-04-07T09:46:26.229481Z","iopub.status.idle":"2023-04-07T09:46:28.656067Z","shell.execute_reply.started":"2023-04-07T09:46:26.229413Z","shell.execute_reply":"2023-04-07T09:46:28.655319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Adding Augmentation\nHere we use augmentation to get more data into the model","metadata":{"_uuid":"610583c12f54b1cefa7df41cc7bd233907abdd88"}},{"cell_type":"code","source":"from keras.preprocessing.image import ImageDataGenerator\ndg_args = dict(featurewise_center = False, \n                  samplewise_center = False,\n                  rotation_range = 5, \n                  width_shift_range = 0.05, \n                  height_shift_range = 0.05, \n                  shear_range = 0.01,\n                  zoom_range = [0.8, 1.2],  \n               # anatomically it doesnt make sense, but many images are flipped\n                  horizontal_flip = True,  \n                  vertical_flip = False,\n                  fill_mode = 'nearest',\n               data_format = 'channels_last')\n\nimage_gen = ImageDataGenerator(**dg_args)\n\ndef gen_augmented_pairs(in_vol, in_seg, batch_size = 16):\n    while True:\n        seed = np.random.choice(range(9999))\n        # keep the seeds syncronized otherwise the augmentation to the images is different from the masks\n        g_vol = image_gen.flow(in_vol, batch_size = batch_size, seed = seed)\n        g_seg = image_gen.flow(in_seg, batch_size = batch_size, seed = seed)\n        for i_vol, i_seg in zip(g_vol, g_seg):\n            yield i_vol, i_seg","metadata":{"_uuid":"9c4261b0867e499cd47a2d7b2ac3659518c39a8e","execution":{"iopub.status.busy":"2023-04-07T09:46:28.657271Z","iopub.execute_input":"2023-04-07T09:46:28.657741Z","iopub.status.idle":"2023-04-07T09:46:28.679887Z","shell.execute_reply.started":"2023-04-07T09:46:28.657685Z","shell.execute_reply":"2023-04-07T09:46:28.679145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_gen = gen_augmented_pairs(train_vol, train_seg, batch_size = 16)\ntest_gen = gen_augmented_pairs(test_vol, test_seg, batch_size = 16)\ntrain_X, train_Y = next(train_gen)\ntest_X, test_Y = next(test_gen)\nprint(train_X.shape, train_Y.shape)\nprint(test_X.shape, test_Y.shape)","metadata":{"_uuid":"287ec6c8593ab1cb647f3a7467b05d930343c3dd","execution":{"iopub.status.busy":"2023-04-07T09:46:28.681077Z","iopub.execute_input":"2023-04-07T09:46:28.681395Z","iopub.status.idle":"2023-04-07T09:46:29.798945Z","shell.execute_reply.started":"2023-04-07T09:46:28.681316Z","shell.execute_reply":"2023-04-07T09:46:29.798228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training Data","metadata":{"_uuid":"968f50d1a1f8a34b928e8b16480a017508081604"}},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (20, 10))\nax1.imshow(montage2d(train_X[:, :, :, 0]), cmap = 'bone')\nax1.set_title('CXR Image')\nax2.imshow(montage2d(train_Y[:, :, :, 0]), cmap = 'bone')\nax2.set_title('Seg Image')","metadata":{"_uuid":"e5338cbfd43643b19635fde6c8903cfba9551377","execution":{"iopub.status.busy":"2023-04-07T09:46:29.800252Z","iopub.execute_input":"2023-04-07T09:46:29.800585Z","iopub.status.idle":"2023-04-07T09:46:30.383250Z","shell.execute_reply.started":"2023-04-07T09:46:29.800536Z","shell.execute_reply":"2023-04-07T09:46:30.382617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Validation Data","metadata":{"_uuid":"c2c1105ce573b856fb4d32bd0b5c29ec44a4f7f3"}},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (20, 10))\nax1.imshow(montage2d(test_X[:, :, :, 0]), cmap = 'bone')\nax1.set_title('CXR Image')\nax2.imshow(montage2d(test_Y[:, :, :, 0]), cmap = 'bone')\nax2.set_title('Seg Image')","metadata":{"_uuid":"802142fdb126d093a37d79a2963da2eabcd003c0","execution":{"iopub.status.busy":"2023-04-07T09:46:30.384475Z","iopub.execute_input":"2023-04-07T09:46:30.384985Z","iopub.status.idle":"2023-04-07T09:46:30.956977Z","shell.execute_reply.started":"2023-04-07T09:46:30.384920Z","shell.execute_reply":"2023-04-07T09:46:30.956075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show Untrained Results\nHere we show random untrained results","metadata":{"_uuid":"a04ed372bf65e2b30f5e40250821586cdcc4b17b"}},{"cell_type":"code","source":"from skimage.segmentation import mark_boundaries\nfrom skimage.color import label2rgb\ntry:\n    from skimage.util.montage import montage2d\nexcept:\n    from skimage.util import montage2d\ndef add_boundary(in_img, in_seg, cmap = 'bone', norm = True, add_labels = True):\n    if norm:\n        n_img = (1.0*in_img-in_img.min())/(1.1*(in_img.max()-in_img.min()))\n    else:\n        n_img = in_img\n    rgb_img = plt.cm.get_cmap(cmap)(n_img)[:, :, :3]\n    if add_labels:\n        return label2rgb(image = rgb_img, label = in_seg.astype(int), bg_label = 0)\n    else:\n        return mark_boundaries(image = rgb_img, label_img = in_seg.astype(int), color = (0, 1, 0), mode = 'thick')\ndef show_full_st(in_img, in_seg, gt_seg):\n    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (30, 10))\n    out_mtg = add_boundary(montage2d(in_img[:, :, :, 0]), \n                           montage2d(gt_seg[:, :, :, 0]>0.5))\n    ax1.imshow(out_mtg)\n    ax1.set_title('Ground Truth')\n    out_mtg = add_boundary(montage2d(in_img[:, :, :, 0]), \n                           montage2d(in_seg[:, :, :, 0]>0.5))\n    ax2.imshow(out_mtg)\n    ax2.set_title('Prediction')\n    out_mtg = montage2d(in_seg[:, :, :, 0]-gt_seg[:, :, :, 0])\n    ax3.imshow(out_mtg, cmap='RdBu', vmin=-1, vmax=1)\n    ax3.set_title('Difference')\ndef show_examples(n=1, with_roi = True):\n    roi_func = lambda x: x[:, \n                               OUT_DIM[0]//2-32:OUT_DIM[0]//2+32,\n                               OUT_DIM[1]//2-64:OUT_DIM[1]//2,\n                               :\n                              ]\n    for (test_X, test_Y), _ in zip(test_gen, range(n)):\n        seg_Y = net.predict(test_X)\n        show_full_st(test_X, seg_Y, test_Y)\n        show_full_st(roi_func(test_X), roi_func(seg_Y), roi_func(test_Y))\n\nshow_examples(1)","metadata":{"_uuid":"57d2e13f2d09ee0b7e1b07a8fa571963916ce48e","execution":{"iopub.status.busy":"2023-04-07T09:46:30.958272Z","iopub.execute_input":"2023-04-07T09:46:30.958661Z","iopub.status.idle":"2023-04-07T09:46:36.880460Z","shell.execute_reply.started":"2023-04-07T09:46:30.958530Z","shell.execute_reply":"2023-04-07T09:46:36.879644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras.callbacks import ModelCheckpoint, LearningRateScheduler, EarlyStopping, ReduceLROnPlateau\nweight_path=\"{}_weights.best.hdf5\".format('cxr_reg')\n\ncheckpoint = ModelCheckpoint(weight_path, monitor='val_loss', verbose=1, \n                             save_best_only=True, mode='min', save_weights_only = True)\n\nreduceLROnPlat = ReduceLROnPlateau(monitor='val_loss', factor=0.5, \n                                   patience=3, \n                                   verbose=1, mode='min', epsilon=0.0001, cooldown=2, min_lr=1e-6)\nearly = EarlyStopping(monitor=\"val_loss\", \n                      mode=\"min\", \n                      patience=15) # probably needs to be more patient, but kaggle time is limited\ncallbacks_list = [checkpoint, early, reduceLROnPlat]","metadata":{"_uuid":"7547559409daafaca582d8402215aee891df5807","execution":{"iopub.status.busy":"2023-04-07T09:46:36.881934Z","iopub.execute_input":"2023-04-07T09:46:36.882518Z","iopub.status.idle":"2023-04-07T09:46:36.910308Z","shell.execute_reply.started":"2023-04-07T09:46:36.882455Z","shell.execute_reply":"2023-04-07T09:46:36.909483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import clear_output\nloss_history = net.fit_generator(train_gen, \n                  steps_per_epoch=len(train_vol)//train_X.shape[0],\n                  epochs = 25,\n                  validation_data = (test_vol, test_seg),\n                  callbacks=callbacks_list\n                 )\nclear_output()","metadata":{"_uuid":"55778a4baa684c73754dfe4073148b3379357306","execution":{"iopub.status.busy":"2023-04-07T09:46:36.912022Z","iopub.execute_input":"2023-04-07T09:46:36.912720Z","iopub.status.idle":"2023-04-07T09:48:13.448197Z","shell.execute_reply.started":"2023-04-07T09:46:36.912654Z","shell.execute_reply":"2023-04-07T09:48:13.447468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"net.load_weights(weight_path)\nnet.save('full_model.h5')","metadata":{"_uuid":"b7d9ea391d4126d08bc4163bca5fb26040fd07fd","execution":{"iopub.status.busy":"2023-04-07T09:48:13.449498Z","iopub.execute_input":"2023-04-07T09:48:13.449774Z","iopub.status.idle":"2023-04-07T09:48:13.687815Z","shell.execute_reply.started":"2023-04-07T09:48:13.449728Z","shell.execute_reply":"2023-04-07T09:48:13.687146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (10, 5))\nax1.plot(loss_history.history['loss'], '-', label = 'Loss')\nax1.plot(loss_history.history['val_loss'], '-', label = 'Validation Loss')\nax1.legend()\n\nax2.plot(100*np.array(loss_history.history['binary_accuracy']), '-', \n         label = 'Accuracy')\nax2.plot(100*np.array(loss_history.history['val_binary_accuracy']), '-',\n         label = 'Validation Accuracy')\nax2.legend()","metadata":{"_uuid":"1fe3a962ee2944a0852262b9def33fc06e0f83f3","execution":{"iopub.status.busy":"2023-04-07T09:48:13.688888Z","iopub.execute_input":"2023-04-07T09:48:13.689146Z","iopub.status.idle":"2023-04-07T09:48:13.907017Z","shell.execute_reply.started":"2023-04-07T09:48:13.689103Z","shell.execute_reply":"2023-04-07T09:48:13.906225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Show results on the training data","metadata":{"_uuid":"b9037c4427194cb15c2b6d2a1c2811480736e740"}},{"cell_type":"code","source":"show_examples(2)","metadata":{"_uuid":"f196fd5dc39ec12387cd6a0baa36ebe3a2dc090d","execution":{"iopub.status.busy":"2023-04-07T09:48:13.908161Z","iopub.execute_input":"2023-04-07T09:48:13.908630Z","iopub.status.idle":"2023-04-07T09:48:24.137276Z","shell.execute_reply.started":"2023-04-07T09:48:13.908574Z","shell.execute_reply":"2023-04-07T09:48:24.136466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Apply to RSNA Data\nHere we load the RSNA data and apply the model to all of the images","metadata":{"_uuid":"64cfe74ab8c82d225077ae86f63154ac8a663a34"}},{"cell_type":"code","source":"import pydicom\nfrom glob import glob\nbase_rsna_dir = os.path.join('..', 'input', 'rsna-pneumonia-detection-challenge')\ntest_mean, test_std = test_X.mean(), test_X.std()\ndef read_dicom_as_float(in_path):\n    out_mat = pydicom.read_file(in_path).pixel_array\n    norm_mat = (out_mat-1.0*np.mean(out_mat))/np.std(out_mat)\n    # make the RSNA distribution look like the training distribution\n    norm_mat = norm_mat*test_std+test_mean\n    return np.expand_dims(norm_mat, -1).astype(np.float32)\nall_rsna_df = pd.DataFrame({'path': glob(os.path.join(base_rsna_dir, \n                                                      'stage_*_images', '*.dcm'))})\nall_rsna_df.sample(3)","metadata":{"_uuid":"d4624ab6087786b30e4602cf7d39db1a319e42bb","execution":{"iopub.status.busy":"2023-04-07T09:48:24.138601Z","iopub.execute_input":"2023-04-07T09:48:24.139056Z","iopub.status.idle":"2023-04-07T09:48:26.115926Z","shell.execute_reply.started":"2023-04-07T09:48:24.138909Z","shell.execute_reply":"2023-04-07T09:48:26.115296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras import layers\nin_shape = read_dicom_as_float(all_rsna_df.iloc[0,0]).shape\nin_img = layers.Input(in_shape, name='DICOMInput')\nscale_factor = (2,2)\nds_dicom = layers.AvgPool2D(scale_factor)(in_img)\nunet_out = net(ds_dicom)\nus_out = layers.UpSampling2D(scale_factor)(unet_out)\nunet_big = Model(inputs=[in_img], outputs=[us_out])\nunet_big.save('big_model.h5')\nunet_big.summary()","metadata":{"_uuid":"9c3ea0258e44f72b0b150a422df0cb175b85627e","execution":{"iopub.status.busy":"2023-04-07T09:48:26.117093Z","iopub.execute_input":"2023-04-07T09:48:26.117410Z","iopub.status.idle":"2023-04-07T09:48:26.448103Z","shell.execute_reply.started":"2023-04-07T09:48:26.117314Z","shell.execute_reply":"2023-04-07T09:48:26.447463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, m_axs = plt.subplots(2, 3, figsize = (10, 8))\nfor c_ax, (_, c_row) in zip(m_axs.flatten(), \n                            all_rsna_df.sample(6).iterrows()):\n    c_img = read_dicom_as_float(c_row['path'])\n    c_seg = unet_big.predict(np.expand_dims(c_img, 0))[0]\n    c_ax.imshow(add_boundary(c_img[:, :, 0], c_seg[:, :, 0]>0.5))","metadata":{"_uuid":"b22f15957262db7fbe32c54b157c6ad78b53eab8","execution":{"iopub.status.busy":"2023-04-07T09:48:26.449411Z","iopub.execute_input":"2023-04-07T09:48:26.449701Z","iopub.status.idle":"2023-04-07T09:48:29.319036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import zipfile as zf\n# from io import BytesIO\n# from PIL import Image\n# batch_size = 12\n# with zf.ZipFile('masks.zip', 'w') as f:\n#     for i, c_rows in tqdm(all_rsna_df.groupby(lambda x: x//batch_size)):\n#         cur_x = np.stack(c_rows['path'].map(read_dicom_as_float), 0)\n#         cur_pred = unet_big.predict(cur_x)>0.5\n#         for out_img, (_, c_row) in zip(cur_pred[:, :, :, 0], c_rows.iterrows()):\n#             arc_name = os.path.relpath(c_row['path'], base_rsna_dir)\n#             arc_name, _ = os.path.splitext(arc_name)\n#             out_pil_obj = Image.fromarray((255*out_img).astype(np.uint8))\n#             out_obj = BytesIO()\n#             out_pil_obj.save(out_obj, format='png')\n#             out_obj.seek(0)\n#             f.writestr('{}.png'.format(arc_name), out_obj.read(), zf.ZIP_STORED)","metadata":{"_uuid":"20d7195defa715ef32055ad5e0c1b868bbce17c6","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"!ls -lh *.zip","metadata":{"_uuid":"3ce40a8d228ffcc2624c402c87029699f89e742a","trusted":true},"execution_count":null,"outputs":[]}]}