{"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":"# HuBMAP - Efficient Sampling Baseline (deepflash2, pytorch, fastai) [sub]\n\n> Submission kernel for model trained with efficient region based sampling. \n\n- Train Notebook: https://www.kaggle.com/matjes/hubmap-efficient-sampling-deepflash2-train\n- Sampling Notebook: https://www.kaggle.com/matjes/hubmap-labels-pdf-0-5-0-25-0-01\n\nRequires deepflash2 (git version), zarr, and segmentation-models-pytorch\n\n\n## Overview\n\n1. Installation and package loading\n2. Helper functions and patches\n3. Configuration\n4. Prediction\n5. Submission\n\n### Versions\n- V12: Minor changes in deepflash2 API to support albumentations (changes `apply`in `DeformationField` slightly, see patch below)\n- V13: Adding prediction threshold 0.3","metadata":{"id":"OcsetTMwKXqC"}},{"cell_type":"markdown","source":"### Installation and package loading","metadata":{}},{"cell_type":"code","source":"# Install deepflash2 and dependencies\n! 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\nimport numpy as np\nimport pandas as pd\nimport os\nimport glob\nimport gc\n\nimport rasterio\nfrom rasterio.windows import Window\n\nimport pathlib\nfrom tqdm.notebook import tqdm\nimport cv2\nfrom keras import backend as K\n\nimport tensorflow as tf\nimport efficientnet as efn\nimport efficientnet.tfkeras\n\nimport sys\nsys.path.append(\"../input/zarrkaggleinstall\")\nsys.path.append(\"../input/segmentation-models-pytorch-install\")\n!pip install -q --no-deps ../input/deepflash2-lfs\nimport cv2, torch, zarr, tifffile, pandas as pd, gc\nfrom fastai.vision.all import *\nfrom deepflash2.all import *\nimport segmentation_models_pytorch as smp\nimport os\nimport pathlib\nfrom tqdm.notebook import tqdm\nimport cv2\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport cv2\nimport glob\nimport gc\nimport rasterio\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Helper functions and patches","metadata":{}},{"cell_type":"code","source":"#https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\n#with transposed mask\ndef rle_encode_less_memory(img):\n    #the image should be transposed\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)\n\ndef load_model_weights(model, file, strict=True):\n    state = torch.load(file, map_location='cpu')\n    stats = state['stats']\n    model_state = state['model']\n    model.load_state_dict(model_state, strict=strict)\n    return model, stats","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Patches for deepflash2 classes, see https://fastcore.fast.ai/basics.html#patch","metadata":{}},{"cell_type":"code","source":"# https://matjesg.github.io/deepflash2/data.html#BaseDataset\n# Handling of different input shapes\n@patch\ndef read_img(self:BaseDataset, *args, **kwargs):\n    \n    data = rasterio.open(args[0], transform = identity, num_threads='all_cpus')\n    if data.count != 3:\n        subdatasets = data.subdatasets\n        layers = []\n        if len(subdatasets) > 0:\n            for i, subdataset in enumerate(subdatasets, 0):\n                layers.append(rasterio.open(subdataset))\n                \n            image = np.zeros((data.shape[0],data.shape[1],3),np.uint8)\n            for i,layer in enumerate(layers):\n                image[:,:,i] =layer.read(1)\n    \n    else:\n        image = data.read([1, 2, 3])\n        image = np.moveaxis(image, 0, -1)\n    \n    return image\n\n# https://matjesg.github.io/deepflash2/data.html#DeformationField\n# Adding normalization (divide by 255)\n@patch\ndef apply(self:DeformationField, data, offset=(0, 0), pad=(0, 0), order=1):\n    \"Apply deformation field to image using interpolation\"\n    outshape = tuple(int(s - p) for (s, p) in zip(self.shape, pad))\n    coords = [np.squeeze(d).astype('float32').reshape(*outshape) for d in self.get(offset, pad)]\n    # Get slices to avoid loading all data (.zarr files)\n    sl = []\n    for i in range(len(coords)):\n        cmin, cmax = int(coords[i].min()), int(coords[i].max())\n        dmax = data.shape[i]\n        if cmin<0: \n            cmax = max(-cmin, cmax)\n            cmin = 0 \n        elif cmax>dmax:\n            cmin = min(cmin, 2*dmax-cmax)\n            cmax = dmax\n            coords[i] -= cmin\n        else: coords[i] -= cmin\n        sl.append(slice(cmin, cmax))    \n    if len(data.shape) == len(self.shape) + 1:\n        \n        ## Channel order change in V12\n        tile = np.empty((*outshape, data.shape[-1]))\n        for c in range(data.shape[-1]):\n            # Adding divide\n            tile[..., c] = cv2.remap(data[sl[0],sl[1], c]/255, coords[1],coords[0], interpolation=order, borderMode=cv2.BORDER_REFLECT)\n    else:\n        tile = cv2.remap(data[sl[0], sl[1]], coords[1], coords[0], interpolation=order, borderMode=cv2.BORDER_REFLECT)\n    return tile","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Configuration","metadata":{}},{"cell_type":"code","source":"class CONFIG():\n    \n    # data paths\n    data_path = Path('../input/hubmap-kidney-segmentation')\n    model_file = '../input/hubmap-deepflash-efficientnetb4/'\n    \n    # deepflash2 dataset (https://matjesg.github.io/deepflash2/data.html#TileDataset)\n    scale = 3 # zoom facor (zoom out)\n    tile_shape = (512, 512)\n    padding = (100,100) # Border overlap for prediction\n    pred_tta = True\n\n    # pytorch model (https://github.com/qubvel/segmentation_models.pytorch)\n    encoder_name = \"efficientnet-b4\"\n    encoder_weights = None\n    in_channels = 3\n    classes = 2\n    \n    # dataloader \n    batch_size = 1\n    \n    # prediction threshold\n    threshold = 0.30\n    \ncfg = CONFIG()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import yaml\nimport pprint\n    \nTHRESHOLD = 0.3 # preds > THRESHOLD\nVOTERS = 0.5\nWINDOW = 1024\nMIN_OVERLAP = 300\nNEW_SIZE = 256\nSUM_PRED = 128\nCHECK=False","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import json\nwith open('../input/hubmap-efficientnet-b4/' + 'metrics.json') as json_file:\n    M = json.load(json_file)\nprint('Model run datetime: '+M['datetime'])\nprint('OOF val_dice_coe: ' + str(M['oof_dice_coe']))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def rle_encode_less_memory_tf(img):\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 make_grid_tf(shape, window=1024, min_overlap=32):\n    \"\"\"\n        Return Array of size (N,4), where N - number of tiles,\n        2nd axis represente slices: 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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models_path_list = [\n    '../input/hubmap-efficientnet-b4/model-fold-2.h5',\n    '../input/hubmap-efficientnet-b6-pseudo/model-fold-0.h5',\n    '../input/hubmap-efficientnetb7-pseudo-labelled/model-fold-1.h5'\n]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fold_models = []\nfor fold_model_path in models_path_list:\n    fold_models.append(tf.keras.models.load_model(fold_model_path,compile = False))\nprint(len(fold_models))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"deepflash_models=[]\n\n#for fold in [1, 2]:\n        \n###################### efficient net b4 #########################\n'''\nmodel = smp.Unet(encoder_name='efficientnet-b6', \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel, stats = load_model_weights(model, '../input/hubmap-deepflash-efficientnet-b6-local/'+f'unet_efficientnet-b6_1.pth')\nbatch_tfms = [Normalize.from_stats(*stats)]\n\nif torch.cuda.is_available():  model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()\n\nmodel = smp.Unet(encoder_name='efficientnet-b4', \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel, stats = load_model_weights(model, '../input/hubmap-single-fold-models-b1b5/'+f'unet_efficientnet-b4_1.pth')\nbatch_tfms = [Normalize.from_stats(*stats)]\nif torch.cuda.is_available():  model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()\n\nmodel = smp.Unet(encoder_name='efficientnet-b5', \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel, stats = load_model_weights(model, '../input/hubmap-single-fold-models-b1b5/'+f'unet_efficientnet-b5_1.pth')\nbatch_tfms = [Normalize.from_stats(*stats)]\nif torch.cuda.is_available():  model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()\n        \n'''\n\nmodel = smp.Unet(encoder_name='efficientnet-b7', \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel, stats = load_model_weights(model, '../input/hubmap-deepflash-effiicentnetb7/'+f'unet_efficientnet-b7.pth')\nbatch_tfms = [Normalize.from_stats(*stats)]\nif torch.cuda.is_available():  model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()\n    \nmodel = smp.Unet(encoder_name='efficientnet-b6', \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel, stats = load_model_weights(model, '../input/hubmap-deepflash-efficientnetb6/'+f'unet_efficientnet-b6.pth')\nbatch_tfms = [Normalize.from_stats(*stats)]\nif torch.cuda.is_available(): model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(deepflash_models))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Prediction","metadata":{}},{"cell_type":"code","source":"import pathlib\nimport glob\nfrom tqdm.notebook import tqdm\np = pathlib.Path('../input/hubmap-kidney-segmentation')\nids=[]\n\nfor index, filename in tqdm(enumerate(p.glob('test/*.tiff')), \n                        total = len(list(p.glob('test/*.tiff')))):\n    ids.append(filename.stem)\n\ndf_sample=pd.DataFrame()\ndf_sample['id']=ids\ndf_sample['predicted']=np.nan\ndf_sample = df_sample.set_index('id')\n\nif df_sample.shape[0]==5:\n    df_sample = df_sample.iloc[:2, :]\nelse:\n    df_sample=df_sample","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from numba import cuda\nnames,preds = [],[]\n\nweights=[1, 2]\n    \nfor idx, _ in df_sample.iterrows():\n    print(f'###### File {idx} ######')\n    \n    # RESTRICT TENSORFLOW TO 2GB OF GPU RAM\n    # SO THAT WE HAVE 14GB RAM FOR RAPIDS\n    LIMIT = 6.0\n    gpus = tf.config.experimental.list_physical_devices('GPU')\n    if gpus:\n        try:\n            tf.config.experimental.set_virtual_device_configuration(\n                gpus[0],\n                [tf.config.experimental.VirtualDeviceConfiguration(memory_limit=1024*LIMIT)])\n            logical_gpus = tf.config.experimental.list_logical_devices('GPU')\n            #print(len(gpus), \"Physical GPUs,\", len(logical_gpus), \"Logical GPUs\")\n        except RuntimeError as e:\n            print(e)\n    print('We will restrict TensorFlow to max %iGB GPU RAM'%LIMIT)\n    print('then RAPIDS can use %iGB GPU RAM'%(16-LIMIT))\n\n    f = cfg.data_path/'test'/f'{idx}.tiff'\n    \n    # Model\n    pred=None\n    \n    dataset_tf = rasterio.open(f, transform = identity)\n    pred_tf = np.zeros(dataset_tf.shape, dtype=np.uint8)  \n    \n    slices_tf = make_grid_tf(dataset_tf.shape, window=WINDOW, min_overlap=MIN_OVERLAP)\n\n    if dataset_tf.count != 3:\n        print('Image file with subdatasets as channels')\n        layers_tf = [rasterio.open(subd) for subd in dataset_tf.subdatasets]\n\n    print(f'Dataset Shape: {dataset_tf.shape}')\n    \n    EMPTY = np.zeros((NEW_SIZE, NEW_SIZE))\n    \n    for (x1,x2,y1,y2) in tqdm(slices_tf):\n            if dataset_tf.count == 3:\n                image = dataset_tf.read([1,2,3],\n                            window=Window.from_slices((x1,x2),(y1,y2)))\n                image = np.moveaxis(image, 0, -1)\n            else:\n                image = np.zeros((WINDOW, WINDOW, 3), dtype=np.uint8)\n                for fl in range(3):\n                    image[:,:,fl] = layers_tf[fl].read(window=Window.from_slices((x1,x2),(y1,y2)))\n\n\n            hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)\n            h,s,v = cv2.split(hsv)\n            s_th = 40\n            p_th = 1000*(1024//256)**2\n\n            if (s>s_th).sum() <= p_th or image.sum() <= p_th :\n                pred_temp = EMPTY\n            else:\n\n                image = cv2.resize(image, (NEW_SIZE, NEW_SIZE),interpolation = cv2.INTER_AREA)\n                image = cv2.cvtColor(image, cv2.COLOR_RGB2BGR)\n                image = np.expand_dims(image, 0)\n                #image = tf.cast(image, tf.float32)/255.0\n\n                pred_temp = None\n\n                for fold_model in fold_models:\n                    if pred_temp is None:\n                        pred_temp = np.squeeze(fold_model.predict(image))\n                    else:\n                        pred_temp += np.squeeze(fold_model.predict(image))\n\n                pred_temp = pred_temp/len(fold_models)\n\n\n            pred_temp = cv2.resize(pred_temp, (WINDOW, WINDOW))\n            pred_tf[x1:x2,y1:y2] +=(pred_temp > THRESHOLD).astype(np.uint8) \n    \n    pred_tf = (pred_tf >= VOTERS).astype(np.uint8)\n    np.save('pred_tf.npy', pred_tf)\n    \n    del pred_tf, EMPTY, slices_tf, dataset_tf, fold_model, pred_temp, image, hsv, h, s, v\n    gc.collect()\n    K.clear_session()\n    \n    \n   \n    # Create deepflash2 dataset (including tiling and file conversion)\n    ds = TileDataset([f], scale=cfg.scale, tile_shape=cfg.tile_shape, padding=cfg.padding)\n    shape = ds.data[f.name].shape\n    print('Shape:', shape)\n\n    \n    # Create fastai dataloader and learner\n    dls = DataLoaders.from_dsets(ds, batch_size=cfg.batch_size, after_batch=batch_tfms, shuffle=False, drop_last=False)\n    \n    for model, weight in zip(deepflash_models, weights):\n        print(weight)\n        if torch.cuda.is_available(): dls.cuda(), model.cuda()\n    \n        learn = Learner(dls, model, loss_func='')\n\n        print('Prediction')\n        if weight==2:\n            res = learn.predict_tiles(dl=dls.train, path='/kaggle/temp/', use_tta=True, uncertainty_estimates=False)\n        else:\n            res = learn.predict_tiles(dl=dls.train, path='/kaggle/temp/', use_tta=True, uncertainty_estimates=False)\n        if pred is None:\n            pred = weight*res[0][f.name][..., 1]/np.sum(weights)\n        else:\n            pred += weight*res[0][f.name][..., 1]/np.sum(weights)\n\n        del res, learn, model\n        gc.collect()\n    \n    del ds, dls\n    gc.collect()\n    torch.cuda.empty_cache()\n    \n    # Load mask from softmax prediction > threshold\n    \n    msk = (pred>cfg.threshold).astype(np.uint8)\n \n    print('Rezising')\n    msk = cv2.resize(msk, (shape[1], shape[0]))\n    pred_tf = np.load('../working/pred_tf.npy')\n    msk = (msk+pred_tf)\n    del pred_tf\n    gc.collect()\n    msk = (msk>=1).astype(np.uint8)\n    rle = rle_encode_less_memory(msk)\n    \n    names.append(idx)\n    preds.append(rle)\n    \n    del msk, rle\n    gc.collect()\n    \n\n    # Overwrite store (reduce disk usage)\n    _ = [shutil.rmtree(p, ignore_errors=True) for p in Path('/kaggle/temp/').iterdir()]\n    _ = [shutil.rmtree(p, ignore_errors=True) for p in Path('/tmp/').iterdir() if p.name.startswith('zarr')]\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Submission","metadata":{}},{"cell_type":"code","source":"df = pd.DataFrame({'id':names,'predicted':preds}).set_index('id')\ndf_sample.loc[df.index.values] = df.values  \ndf_sample.to_csv('submission.csv')\ndisplay(df_sample)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}