{"cells":[{"metadata":{},"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 [[data] HuBMAP Image 2 TFRecords 128,256,512,1024](https://www.kaggle.com/mistag/data-hubmap-image-2-tfrecords-128-256-512-1024). The TFRecords are also available in a [dataset](https://www.kaggle.com/mistag/hubmap-tfrecords) (needed for TPU training).   \n\nTraining of the model is done in [this notebook](https://www.kaggle.com/mistag/train-fpn-efficientnetb2) (not public yet!), with validation on one image and training on the other 7 in a K-fold cross-validation scheme.\n"},{"metadata":{},"cell_type":"markdown","source":"First install a few libraries needed with Segmentation Models."},{"metadata":{"trusted":true},"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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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\nfrom skimage import io\nimport json","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## File functions\nInput images are really big, and to save memory the files are mapped to disk with numpy.memmap(). A bit slow, but frees a lot of memory."},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_tif_file(fname):\n    img = io.imread(fname)\n    img = np.squeeze(img)\n    if img.shape[0] == 3: # swap axes as required\n        img = img.swapaxes(0,1)\n        img = img.swapaxes(1,2)\n    return img\n\n# map image to file(s)\ndef map_img2file(fname):\n    img = read_tif_file(fname)\n    dims = np.array(img.shape)\n    ch = 1 if len(dims) == 2 else dims[2]\n    for i in range(ch):\n        f = np.memmap('img{}.dat'.format(i), dtype=np.uint8, mode='w+', shape=(dims[0], dims[1]))\n        f[:] = img[:,:,i] if ch > 1 else img[:,:]\n        del f\n    return dims\n\n# read part of image from file\ndef get_patch_from_file(dims, pos, psize):\n    ch = 1 if len(dims) == 2 else dims[2]\n    patch = np.zeros([psize[0], psize[1]], dtype=np.uint8) if ch == 1 else np.zeros([psize[0], psize[1], ch], dtype=np.uint8)\n    for i in range(ch):\n        f = np.memmap('img{}.dat'.format(i), dtype=np.uint8, mode='r', shape=(dims[0], dims[1]))\n        p = f[pos[0]:pos[0]+psize[0], pos[1]:pos[1]+psize[1]]\n        crop = p.shape\n        if ch == 1:\n            patch[0:p.shape[0], 0:p.shape[1]] = p\n        else:\n            patch[0:p.shape[0], 0:p.shape[1],i] = p\n        del f\n    return patch, crop","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using the RLE-encoder from [this notebook](https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub):"},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## TTA\nTest-time augmentation functions."},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{},"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. The mask of image afa5e8098 is shifted according to [this discussion](https://www.kaggle.com/c/hubmap-kidney-segmentation/discussion/207517). 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":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"MEANING_OF_LIFE = 42\nMEANING_OF_LIFE_REV = int(str(MEANING_OF_LIFE)[::-1])\n\nPATH = '../input/hubmap-kidney-segmentation/test/'\nfilelist = glob.glob(PATH+'*.tiff')\nSUB_FILE = './submission.csv'\nwith open(SUB_FILE, 'w') as f:\n    f.write(\"id,predicted\\n\")\n    \nMODELS = ['../input/train-fpn-efficientnetb2/FPN+ENetB4-2']\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    print(fid)\n    dims = map_img2file(file)\n    pmask = np.zeros(dims[:2], dtype=np.uint8)\n    x_shft, y_shft = 0,0\n    if fid == 'afa5e8098': # mask correction\n        x_shft, y_shft = MEANING_OF_LIFE, MEANING_OF_LIFE_REV\n    for modl in range(len(MODELS)):\n        print(MODELS[modl])\n        # load pre-trained model\n        mname = MODELS[modl]\n        with open(mname+'.json', 'r') as m:\n            lm = m.read()\n            model = model_from_json(lm)\n        model.load_weights(mname+'.h5')\n        # process image\n        for x in range((dims[0]-OVERLAP-x_shft)//STEP + min(1,(dims[0]-OVERLAP-x_shft) % STEP)):\n            for y in range((dims[1]-OVERLAP-y_shft)//STEP + min(1,(dims[1]-OVERLAP-y_shft) % STEP)):\n                tile, crop = get_patch_from_file(dims, [x*STEP+x_shft, y*STEP+y_shft], [size,size])\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_BGR2HSV))\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 >= MEANING_OF_LIFE_REV - K_SPLITS\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\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# clean up intermediate files\n%rm -f *.dat","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}