{"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":"markdown","source":"# Submission\n\nThis notebook contains the inference code for the HuBMAP competition. For our submission, we used a U-Net model with the Tversky loss function and EfficientNet backbone. The model was pretrained on the ImageNet dataset. In order to avoid OOM, images are split into frames on the fly and processed in small batches. In the end the resulting mask. Our submission implements the following steps:\n1. Extraction of frames from the original image\n2. Downscaling of frames according to the input size of the network\n3. Pixel normalization\n4. Inference\n5. Upscaling to the original size of a frame\n6. Cropping the center of a frame in order to avoid the edge effect (optional)\n7. Conversion of the inferred masks to binary masks using threshold\n8. Reconstruction of the original image\n9. Conversion of binary masks to RLE-encoded strings\n10. Plotting the resulting masks (optional)\n11. Storing RLE-encoded masks on the file system\n\nSetup libraries that are necessary for the submission.","metadata":{}},{"cell_type":"code","source":"!pip install ../input/keras-applications/Keras_Applications-1.0.8/ -f ./ --no-index\n!pip install ../input/image-classifiers/image_classifiers-1.0.0/ -f ./ --no-index\n!pip install ../input/efficientnet-1-0-0/efficientnet-1.0.0/ -f ./ --no-index\n!pip install ../input/segmentation-models/segmentation_models-1.0.1/ -f ./ --no-index\n%env SM_FRAMEWORK=tf.keras","metadata":{"_kg_hide-output":true,"_kg_hide-input":false,"execution":{"iopub.status.busy":"2022-05-04T14:42:32.832657Z","iopub.execute_input":"2022-05-04T14:42:32.833071Z","iopub.status.idle":"2022-05-04T14:43:01.864585Z","shell.execute_reply.started":"2022-05-04T14:42:32.832984Z","shell.execute_reply":"2022-05-04T14:43:01.863369Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import cv2\nimport gc\nimport glob\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport rasterio\nimport segmentation_models as sm\nimport tensorflow as tf\nimport tifffile as tiff\nfrom rasterio.windows import Window\nfrom tensorflow import keras\nfrom tensorflow.keras import backend as K","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-04T14:43:01.866334Z","iopub.execute_input":"2022-05-04T14:43:01.866712Z","iopub.status.idle":"2022-05-04T14:43:07.012044Z","shell.execute_reply.started":"2022-05-04T14:43:01.866671Z","shell.execute_reply":"2022-05-04T14:43:07.011131Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our solution uses the Tversky loss function that needs to be defined manually since it is not part of the tensorflow library. The implementation of the loss function was taken from https://www.kaggle.com/code/bigironsphere/loss-function-library-keras-pytorch. ","metadata":{}},{"cell_type":"code","source":"ALPHA = 0.5\nBETA = 0.5\n\ndef TverskyLoss(targets, inputs, alpha=ALPHA, beta=BETA, smooth=1e-6):\n    '''\n    source: https://www.kaggle.com/code/bigironsphere/loss-function-library-keras-pytorch\n    \n    tversky loss function for keras\n    '''\n    \n    #flatten label and prediction tensors\n    inputs = K.flatten(inputs)\n    targets = K.flatten(targets)\n\n    #True Positives, False Positives & False Negatives\n    TP = K.sum((inputs * targets))\n    FP = K.sum(((1-targets) * inputs))\n    FN = K.sum((targets * (1-inputs)))\n\n    Tversky = (TP + smooth) / (TP + alpha*FP + beta*FN + smooth)  \n\n    return 1 - Tversky","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:07.014050Z","iopub.execute_input":"2022-05-04T14:43:07.014388Z","iopub.status.idle":"2022-05-04T14:43:07.021108Z","shell.execute_reply.started":"2022-05-04T14:43:07.014350Z","shell.execute_reply":"2022-05-04T14:43:07.020297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Set notebook variables.","metadata":{}},{"cell_type":"code","source":"# general notebook variables\nDATASET_PATH = '/kaggle/input/hubmap-kidney-segmentation'\nTEST_DATASET_PATH = os.path.join(DATASET_PATH, 'test')\nTEST_IMG_IDS = [os.path.basename(img_id).split(\".\")[0] for img_id in glob.glob(f'{TEST_DATASET_PATH}/*.tiff')]\nMODEL_PATH = '../input/notest-trainweights-unet-efficentnetb0/unet_efficientnetb0_f01_1.h5'\nSUB_CSV = pd.read_csv(os.path.join(DATASET_PATH, 'sample_submission.csv'), index_col='id')\nREMOVE_EDGES = False # remove edges to avoid the edge effect\nPLOT_RESULT = False # plot resulting masks\n\n# frame variables\nINPUT_SIZE = (512, 512) # network input size\nFRAME_SIZE = 1024 # size of a frame\nFRAME_OVERLAP = 384 # frame overlap size\nFRAME_DISCARD = FRAME_OVERLAP // 2 # size of the frame edges that will be removed if REMOVE_EDGES = True\nMASK_THRESH = 0.4\nBATCH_SIZE = 16\n\n# keras settings\nLR = 5e-4\nOPTIMIZER = keras.optimizers.RMSprop(LR)\nLOSS_FUNC = TverskyLoss\nMETRICS_ARR = [sm.metrics.f1_score, sm.metrics.iou_score]","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:07.022969Z","iopub.execute_input":"2022-05-04T14:43:07.023641Z","iopub.status.idle":"2022-05-04T14:43:07.052078Z","shell.execute_reply.started":"2022-05-04T14:43:07.023564Z","shell.execute_reply":"2022-05-04T14:43:07.051324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Load weights of a trained model and compile.","metadata":{}},{"cell_type":"code","source":"model = keras.models.load_model(MODEL_PATH, compile=False)\nmodel.compile(OPTIMIZER, loss=LOSS_FUNC, metrics=METRICS_ARR) # same params as in train","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:07.053224Z","iopub.execute_input":"2022-05-04T14:43:07.053558Z","iopub.status.idle":"2022-05-04T14:43:11.541802Z","shell.execute_reply.started":"2022-05-04T14:43:07.053525Z","shell.execute_reply":"2022-05-04T14:43:11.540957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Define auxiliary functions. The function `make_grid` (borrowed from https://www.kaggle.com/leighplt/pytorch-fcn-resnet50) creates a grid for the sliding window operation. `rle_encode_less_memory` function (from https://www.kaggle.com/code/bguberfain/memory-aware-rle-encoding) converts the image to an RLE-encoded string. ","metadata":{}},{"cell_type":"code","source":"def make_grid(shape, window=256, min_overlap=32):\n    \"\"\"\n    source: https://www.kaggle.com/leighplt/pytorch-fcn-resnet50\n    \n    function to generate a grid layout for sliding window\n    :param shape: height and width of the image\n    :param window: size of the window\n    :param min_overlap: minimal window overlap\n    :return: array of window coordinates (x1,x2,y1,y2)\n    \"\"\"\n    x, y = shape\n    nx = x // (window - min_overlap) + 1\n    x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)\n    x1[-1] = x - window\n    x2 = (x1 + window).clip(0, x)\n    ny = y // (window - min_overlap) + 1\n    y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)\n    y1[-1] = y - window\n    y2 = (y1 + window).clip(0, y)\n    slices = np.zeros((nx,ny, 4), dtype=np.int64)\n    \n    for i in range(nx):\n        for j in range(ny):\n            slices[i,j] = x1[i], x2[i], y1[j], y2[j]    \n    return slices.reshape(nx*ny,4)\n\ndef rle_encode_less_memory(img):\n    \"\"\"\n    source: https://www.kaggle.com/code/bguberfain/memory-aware-rle-encoding\n    \n    function This method requires first and last pixel to be zero\n    :param img: numpy array, 1 - mask, 0 - background\n    :return: run length as string formated\n    \"\"\"\n    pixels = img.T.flatten()\n    pixels[0] = 0\n    pixels[-1] = 0\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 2\n    runs[1::2] -= runs[::2]\n    return ' '.join(str(x) for x in runs)\n\ndef remove_edges(frame, coordinates, discard):\n    \"\"\"\n    function to crop the center of a frame\n    :param frame: frame to be processed\n    :param coordinates: coordinates of the frame in the original image\n    :param discard: size of margin to be removed\n    :return: cropped frame and its new coordinates\n    \"\"\"\n    x1,x2,y1,y2 = coordinates\n    \n    if x1 != 0:\n        x1 += discard\n        frame = frame[:,discard:]\n    if x2 != shape[1]:\n        x2 -= discard\n        frame = frame[:,:-discard or None]\n    if y1 !=0:\n        y1 += discard\n        frame = frame[discard:,:]\n    if y2 != shape[0]:\n        y2 -= discard\n        frame = frame[:-discard or None,:]\n    return (frame, (x1,x2,y1,y2))\n\ndef plot_result(dataset_path, img_id, mask):\n    \"\"\"\n    function to plot inferred masks\n    :param dataset_path: path to the test dataset\n    :param img_id: name of the image file (without extension)\n    :param mask: inferred mask\n    :return: None\n    \"\"\"\n    full_image = tiff.imread(os.path.join(dataset_path, f'{img_id}.tiff'))\n        \n    # reshape image if necessary\n    if len(full_image.shape) == 5:\n        full_image = full_image.squeeze()\n    if full_image.shape[0] == 3:\n        full_image = full_image.transpose(1, 2, 0)\n\n    # create a copy of the image\n    full_image = full_image.copy()\n\n    scale_factor = 10\n    full_image = cv2.resize(full_image,(full_image.shape[1]//scale_factor,full_image.shape[0]//scale_factor), interpolation = cv2.INTER_CUBIC)\n    mask = cv2.resize(mask,(mask.shape[1]//scale_factor, mask.shape[0]//scale_factor), interpolation = cv2.INTER_CUBIC)\n\n    plt.figure(figsize=(15,15))\n    plt.imshow(full_image)\n    plt.imshow(mask, alpha=0.5, cmap='plasma')\n    plt.show()\n\n    del full_image\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:11.543064Z","iopub.execute_input":"2022-05-04T14:43:11.543383Z","iopub.status.idle":"2022-05-04T14:43:11.561556Z","shell.execute_reply.started":"2022-05-04T14:43:11.543350Z","shell.execute_reply":"2022-05-04T14:43:11.560799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Optionally, our inference code crops the center of every frame in order to avoid the edge effect as described in https://www.kaggle.com/competitions/hubmap-kidney-segmentation/discussion/238013.","metadata":{}},{"cell_type":"code","source":"for img_id in TEST_IMG_IDS:\n    print(img_id)\n    \n    # source: https://www.kaggle.com/code/iafoss/256x256-images\n    # some images have issues with their format\n    # and must be saved correctly before reading with rasterio\n    img = rasterio.open(os.path.join(TEST_DATASET_PATH, f'{img_id}.tiff'), num_threads='all_cpus')\n    if img.count != 3: # the number of raster bands in the dataset\n        subdatasets = img.subdatasets # sequence of subdatasets\n        layers = []\n        if len(subdatasets) > 0:\n            for i, subdataset in enumerate(subdatasets, 0):\n                layers.append(rasterio.open(subdataset))\n                \n    # calculate frame coordinates\n    shape = img.shape\n    frames = make_grid((shape[1], shape[0]), window=FRAME_SIZE, min_overlap=FRAME_OVERLAP)\n        \n    # placeholder matrix for the resulting mask\n    wsi_mask = np.ones(shape, np.uint8)\n        \n    # process frames as small batches\n    for batch_idx in range(0, math.ceil(len(frames) / BATCH_SIZE)):\n                \n        batch = [] # list of image frames\n        batch_coordinates = [] # list of frame coordinates\n        \n        # create batch\n        for frame_idx in range(batch_idx * BATCH_SIZE, min((batch_idx+1) * BATCH_SIZE, len(frames))):            \n            img_frame = np.zeros((FRAME_SIZE, FRAME_SIZE, 3), np.uint8)\n            x1,x2,y1,y2 = frames[frame_idx]\n            \n            # source: https://www.kaggle.com/code/iafoss/256x256-images\n            if img.count == 3:\n                img_frame = np.moveaxis(img.read([1,2,3], window=Window.from_slices((y1,y2),(x1,x2))), 0, -1)\n            else:\n                for i,layer in enumerate(layers):\n                    img_frame[:,:,i] = layer.read(1, window=Window.from_slices((y1,y2),(x1,x2)))\n                                \n            # resize according to the input size of the network\n            img_frame = cv2.resize(img_frame, INPUT_SIZE, interpolation=cv2.INTER_CUBIC)\n            \n            # convert to float values\n            img_frame = tf.cast(img_frame/255.0, tf.float32)\n                        \n            batch.append(img_frame)\n            batch_coordinates.append(frames[frame_idx])\n                \n        # predict batch\n        result = model.predict(np.array(batch))\n        \n        # process batch\n        for res_idx in range(len(result)):\n            # get predicted frame mask\n            res_frame = result[res_idx]\n            res_coordinates = batch_coordinates[res_idx]\n            \n            # resize to original size\n            res_frame = cv2.resize(res_frame, (FRAME_SIZE, FRAME_SIZE), interpolation=cv2.INTER_CUBIC)\n            \n            if REMOVE_EDGES:\n                res_frame,res_coordinates = remove_edges(res_frame, res_coordinates, FRAME_DISCARD)\n            \n            # apply threshold to the predicted mask\n            res_frame = (res_frame > MASK_THRESH).astype(np.uint8)\n            \n            # insert predicted frame into\n            x1,x2,y1,y2 = res_coordinates # get coordinates of the frame\n            wsi_mask[y1:y2,x1:x2] = wsi_mask[y1:y2,x1:x2] * res_frame\n            \n        del res_frame, res_coordinates, result, batch, batch_coordinates, img_frame\n        gc.collect()\n        \n    # save as csv\n    SUB_CSV.loc[img_id,'predicted'] = rle_encode_less_memory(wsi_mask)\n    \n    if PLOT_RESULT:\n        plot_result(TEST_DATASET_PATH, img_id, wsi_mask)\n        \n    del img, frames, wsi_mask\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:11.562924Z","iopub.execute_input":"2022-05-04T14:43:11.563468Z","iopub.status.idle":"2022-05-04T14:43:50.321911Z","shell.execute_reply.started":"2022-05-04T14:43:11.563431Z","shell.execute_reply":"2022-05-04T14:43:50.319912Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Write resulting masks to a CSV file.","metadata":{}},{"cell_type":"code","source":"SUB_CSV.to_csv('submission.csv')\nSUB_CSV","metadata":{"execution":{"iopub.status.busy":"2022-05-04T14:43:50.323128Z","iopub.status.idle":"2022-05-04T14:43:50.323865Z"},"trusted":true},"execution_count":null,"outputs":[]}]}