{"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":"![title](https://camo.githubusercontent.com/51eea85ed59f27be0485cc5774d09b522ea8e77cd3f0753c085cacd18d4a41a0/68747470733a2f2f692e6962622e636f2f4774784753386d2f5365676d656e746174696f6e2d4d6f64656c732d56312d536964652d332d312e706e67)","metadata":{}},{"cell_type":"markdown","source":"# Inference with Keras Segmentation Models Library\nThis notebook will make predictions on the HuBMAP data with a FPN model from the [Segmentation Models library](https://github.com/qubvel/segmentation_models). This library is Keras based and really simple to use. It has four different segmentation models (Unet, Linknet, FPN and PSPNet), and a whopping 25 different pretrained backbones that can be used with each model.  \n\nThe training data has been converted into TFRecords in [HuBMAP Image 2 TFRecords 256,512,1024](https://www.kaggle.com/mistag/hubmap-image-2-tfrecords-256-512-1024). \n\nTraining of the model is done in [this notebook](https://www.kaggle.com/mistag/train-fpn-segmentation-model-no-43e), with validation on three images and training on the other 12 in a K-fold cross-validation scheme.\n\nThe Feature Pyramid Model was used in training:\n![FPN](https://raw.githubusercontent.com/qubvel/segmentation_models/master/images/fpn.png)\n","metadata":{}},{"cell_type":"markdown","source":"First install a few libraries needed with Segmentation Models.","metadata":{}},{"cell_type":"code","source":"!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index -q\n!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index -q","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nos.environ['SM_FRAMEWORK'] = 'tf.keras'\nimport sys\nimport numpy as np\nimport cv2\nimport glob\nimport tensorflow as tf\nfrom tensorflow.keras.layers import *\nfrom tensorflow.keras.models import model_from_json\nfrom tensorflow.keras.utils import CustomObjectScope\nfrom tensorflow.keras import backend as K\nfrom tensorflow.keras.utils import get_custom_objects\nimport efficientnet as efn\nimport efficientnet.tfkeras\nimport rasterio\nfrom rasterio.windows import Window\nimport gdal\nimport json\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get hyperparameters from the training notebook\nwith open('../input/train-fpn-segmentation-model-no-43e/hparams.json') as json_file:\n    hparams = json.load(json_file)\nhparams","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"IMG_SIZE = hparams['IMG_SIZE']\nSCALE_FACTOR = hparams['SCALE_FACTOR']\nK_SPLITS = hparams['K_SPLITS']\nIDNT = rasterio.Affine(1, 0, 0, 0, 1, 0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the RLE-encoder from [this notebook](https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub):","metadata":{}},{"cell_type":"code","source":"##https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\ndef rle_encode_less_memory(img):\n    pixels = img.T.flatten()\n    \n    # This simplified method requires first and last pixel to be zero\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    \n    return ' '.join(str(x) for x in runs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## TTA\nTest-time augmentation functions.","metadata":{}},{"cell_type":"code","source":"def create_TTA_batch(img):\n    if len(img.shape) < 4:\n        img = np.expand_dims(img, 0)\n            \n    batch=np.zeros((img.shape[0]*8,img.shape[1],img.shape[2],img.shape[3]), dtype=np.float32)     \n    for i in range(img.shape[0]):\n        orig = tf.keras.preprocessing.image.img_to_array(img[i,:,:,:])/255. # un-augmented\n        batch[i*8,:,:,:] = orig\n        batch[i*8+1,:,:,:] = np.rot90(orig, axes=(0, 1), k=1)\n        batch[i*8+2,:,:,:] = np.rot90(orig, axes=(0, 1), k=2)\n        batch[i*8+3,:,:,:] = np.rot90(orig, axes=(0, 1), k=3)\n        orig = orig[:, ::-1]\n        batch[i*8+4,:,:,:] = orig\n        batch[i*8+5,:,:,:] = np.rot90(orig, axes=(0, 1), k=1)\n        batch[i*8+6,:,:,:] = np.rot90(orig, axes=(0, 1), k=2)\n        batch[i*8+7,:,:,:] = np.rot90(orig, axes=(0, 1), k=3)\n    return batch\n\ndef mask_TTA(masks):\n    batch=np.zeros((masks.shape[0],masks.shape[1],masks.shape[2],masks.shape[3]), dtype=np.float32)\n    for i in range(masks.shape[0]//8):\n        batch[i*8,:,:,:] = masks[i*8]\n        batch[i*8+1,:,:,:] = np.rot90(masks[i*8+1], axes=(0, 1), k=3)\n        batch[i*8+2,:,:,:] = np.rot90(masks[i*8+2], axes=(0, 1), k=2)\n        batch[i*8+3,:,:,:] = np.rot90(masks[i*8+3], axes=(0, 1), k=1)\n        batch[i*8+4,:,:,:] = masks[i*8+4][:, ::-1]\n        batch[i*8+5,:,:,:] = np.rot90(masks[i*8+5], axes=(0, 1), k=3)[:, ::-1]\n        batch[i*8+6,:,:,:] = np.rot90(masks[i*8+6], axes=(0, 1), k=2)[:, ::-1]\n        batch[i*8+7,:,:,:] = np.rot90(masks[i*8+7], axes=(0, 1), k=1)[:, ::-1]\n    return(batch)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference\nThe images are processed with an overlap of half the patch size, which means that pixels are processed four times. To make an ensemble of several models, just add them to the MODELS list. The mask threshold must be adjusted when adding more models though.","metadata":{}},{"cell_type":"code","source":"MEANING_OF_LIFE = 42\n\nPATH = '../input/hubmap-kidney-segmentation/test/'\nfilelist = glob.glob(PATH+'*.tiff')\nif len(filelist) == 5: # save time durimg commit phase by only processing one image\n    filelist = filelist[:1]\n    COMMIT = True\nelse:\n    COMMIT = False\n\nSUB_FILE = './submission.csv'\nwith open(SUB_FILE, 'w') as f:\n    f.write(\"id,predicted\\n\")\n    \nMODELS = ['../input/train-fpn-segmentation-model-no-43e/FPN-model43e-4']\ns_th = MEANING_OF_LIFE+K_SPLITS # saturation blanking threshold\np_th = IMG_SIZE*IMG_SIZE//32   # pixel count threshold\nsize = int(IMG_SIZE * SCALE_FACTOR)\n\nOVERLAP = size//2\nSTEP = size-OVERLAP\n\nfor file in filelist:\n    fid = file.replace('\\\\','.').replace('/','.').split('.')[-2]\n    img_data = rasterio.open(file, transform=IDNT)\n    if img_data.count != 3: # channels as subdata\n        layers = [rasterio.open(subd) for subd in img_data.subdatasets]\n    dims = [img_data.shape[0], img_data.shape[1]]\n    pmask = np.zeros(dims[:2], dtype=np.uint8)\n    for modl in MODELS:\n        print('image: {}, model: {}'.format(fid, modl))\n        # load pre-trained model\n        with open(modl+'.json', 'r') as m:\n            lm = m.read()\n            model = model_from_json(lm)\n        model.load_weights(modl+'.h5')\n        # process image\n        for x in tqdm(range((dims[0]-OVERLAP)//STEP + min(1,(dims[0]-OVERLAP) % STEP))):\n            for y in range((dims[1]-OVERLAP)//STEP + min(1,(dims[1]-OVERLAP) % STEP)):\n                x1, x2, y1, y2 = x*STEP, x*STEP+size, y*STEP, y*STEP+size\n                crop = [size, size]\n                if x2 > dims[0]:\n                    crop[0] = dims[0] - x*STEP\n                if y2 > dims[1]:\n                    crop[1] = dims[1] - y*STEP\n                \n                if img_data.count == 3: # normal\n                    tile = img_data.read([1, 2, 3], window=Window.from_slices((x1, x2), (y1, y2)))\n                    tile = np.moveaxis(tile, 0, -1)\n                else: # with subdatasets/layers\n                    tile = np.zeros((crop[0], crop[1], 3), dtype=np.uint8)\n                    for fl in range(3):\n                        tile[:, :, fl] = layers[fl].read(window=Window.from_slices((x1, x2), (y1, y2)))\n                if crop != (size,size):\n                    impad = np.zeros((size,size,3), dtype=np.uint8)\n                    impad[:tile.shape[0],:tile.shape[1]] = tile\n                    tile = impad\n                # downscale tile\n                patch = cv2.resize(tile,\n                                   dsize=(IMG_SIZE, IMG_SIZE),\n                                   interpolation = cv2.INTER_AREA)\n                # simple check of saturation if prediction is worthwhile\n                _, s, _ = cv2.split(cv2.cvtColor(patch, cv2.COLOR_RGB2HSV))\n                if (s>s_th).sum() > p_th:\n                    batch = create_TTA_batch(patch)\n                    preds = model.predict(batch)\n                    pred = mask_TTA(preds)\n                    mask = np.rint(np.sum(pred, axis=0))\n                    # upscale tile mask before adding to total mask\n                    pint =cv2.resize(mask.astype(int), dsize=(size, size), interpolation = cv2.INTER_NEAREST)\n                    pmask[x*STEP:x*STEP+crop[0], y*STEP:y*STEP+crop[1]] += pint[0:crop[0], 0:crop[1]].astype(np.uint8)\n    pmask = pmask >= 16\n    with open(SUB_FILE, 'a') as f:\n        f.write(\"{},\".format(fid))\n        f.write(rle_encode_less_memory(pmask))\n        f.write(\"\\n\")","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Reality check: Plot the last pmask.","metadata":{}},{"cell_type":"code","source":"if COMMIT:\n    plt.figure(figsize=(20,20))\n    plt.imshow(pmask);","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}