{"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***\n\n## Overview\n\n1. Installation and package loading\n2. Functions and classes for prediction\n3. Configuration\n4. Prediction\n5. Submission\n\n#### Related Kernels\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\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.4\n- V14: Threshold 0.2 for d488c759a - see discussion https://www.kaggle.com/c/hubmap-kidney-segmentation/discussion/228993 \n- V15: **NEW PREDICTION** \n    - Using overlapping tiles and gaussian weighting from [nnunet](https://www.nature.com/articles/s41592-020-01008-z)/[github](https://github.com/MIC-DKFZ/nnUNet), which will also be part of the upcoming `deepflash2` release\n    - Supporting model ensembles\n    - Fixing submission to private LB using `rasterio` (thanks to @leighplt [kernel](https://www.kaggle.com/leighplt/pytorch-fcn-resnet50) and @iafoss ([kernel](https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub))","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\n\nimport sys\nsys.path.append(\"../input/segmentation-models-pytorch-install\")\n!pip install -q --no-deps ../input/deepflash2-lfs\nimport cv2, torch, gc, rasterio\nimport torch.nn.functional as F\nimport deepflash2.tta as tta\nimport matplotlib.pyplot as plt\nimport pandas as pd, numpy as np\nimport segmentation_models_pytorch as smp\nfrom pathlib import Path\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy .ndimage.filters import gaussian_filter\nfrom tqdm.notebook import tqdm\n\nimport os\nimport glob\nimport gc\n\nimport rasterio\nfrom rasterio.windows import Window\n\nimport pathlib\nfrom tqdm.notebook import tqdm\nimport cv2\n\nimport tensorflow as tf\nimport efficientnet as efn\nimport efficientnet.tfkeras\nimport tensorflow as tf\nfrom tensorflow.keras import backend as K\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Functions and classes for prediction","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\n\n# from https://github.com/MIC-DKFZ/nnUNet/blob/2fade8f32607220f8598544f0d5b5e5fa73768e5/nnunet/network_architecture/neural_network.py#L250\ndef _get_gaussian(patch_size, sigma_scale=1. / 8) -> np.ndarray:\n    tmp = np.zeros(patch_size)\n    center_coords = [i // 2 for i in patch_size]\n    sigmas = [i * sigma_scale for i in patch_size]\n    tmp[tuple(center_coords)] = 1\n    gaussian_importance_map = gaussian_filter(tmp, sigmas, 0, mode='constant', cval=0)\n    gaussian_importance_map = gaussian_importance_map / np.max(gaussian_importance_map) * 1\n    gaussian_importance_map = gaussian_importance_map.astype(np.float32)\n\n    # gaussian_importance_map cannot be 0, otherwise we may end up with nans!\n    gaussian_importance_map[gaussian_importance_map == 0] = np.min(\n        gaussian_importance_map[gaussian_importance_map != 0])\n\n    return gaussian_importance_map","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Some code adapted from https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub\nclass HubmapDataset(Dataset):\n    'HubmapDataset class that does not load the full tiff files.'\n    def __init__(self, file, stats, scale=3, shift=.8, output_shape=(512,512), s_th = 40):\n        \n        self.mean, self.std = stats\n        self.scale = scale\n        self.shift = shift\n        self.output_shape = output_shape\n        self.input_shape = tuple(int(t*scale) for t in self.output_shape)      \n        self.s_th = s_th #saturation blancking threshold\n        self.p_th = 1000*(self.output_shape[0]//256)**2 #threshold for the minimum number of pixels\n\n        identity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n        self.data = rasterio.open(file, transform = identity, num_threads='all_cpus')\n        if self.data.count != 3:\n            subdatasets = self.data.subdatasets\n            self.layers = []\n            if len(subdatasets) > 0:\n                for i, subdataset in enumerate(subdatasets, 0):\n                    self.layers.append(rasterio.open(subdataset))\n            \n        # Tiling\n        self.slices = []\n        self.out_slices = []\n        self.out_data_shape = tuple(int(x//self.scale) for x in self.data.shape)\n        start_points = [o//2 for o in self.output_shape]\n        end_points = [(s - st) for s, st in zip(self.out_data_shape, start_points)]\n        n_points = [int(s//(o*self.shift))+1 for s, o in zip(self.out_data_shape, self.output_shape)]\n        center_points = [np.linspace(st, e, num=n, endpoint=True, dtype=np.int64) for st, e, n in zip(start_points, end_points, n_points)]\n        for cx in center_points[1]:\n            for cy in center_points[0]:\n                # Calculate output slices for whole image\n                slices = tuple(slice(int((c*self.scale - o/2).clip(0, s)), int((c*self.scale + o/2).clip(max=s)))\n                                 for (c, o, s) in zip((cy, cx), self.input_shape, self.data.shape))\n                self.slices.append(slices)\n                \n                out_slices = tuple(slice(int((c - o/2).clip(0, s)), int((c + o/2).clip(max=s)))\n                                 for (c, o, s) in zip((cy, cx), self.output_shape, self.out_data_shape))\n                self.out_slices.append(out_slices)\n                \n\n    def __len__(self):\n        return len(self.slices)\n\n    def __getitem__(self, idx):\n        if torch.is_tensor(idx):\n            idx = idx.tolist()\n            \n        slices = self.slices[idx]\n        if self.data.count == 3: # normal\n            img = self.data.read([1, 2, 3], \n                window=Window.from_slices(*slices)\n            )\n            img = np.moveaxis(img, 0, -1)\n        else: # with subdatasets/layers\n            img = np.zeros((*self.input_shape, 3), dtype=np.uint8)\n            for fl in range(3):\n                img[:, :, fl] = self.layers[fl].read(\n                    window=Window.from_slices(*slices)\n                )\n        \n        if self.scale!=1:\n            img = cv2.resize(img, self.output_shape, interpolation = cv2.INTER_AREA)\n        \n        #check for empty imges\n        hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n        h,s,v = cv2.split(hsv)\n        if (s>self.s_th).sum() <= self.p_th or img.sum() <= self.p_th:\n            # Remove if idx=-1\n            idx = -1\n        \n        img = (img/255.0 - self.mean)/self.std\n        img = img.transpose(2, 0, 1).astype('float32')\n        \n        return torch.from_numpy(img), idx\n    \nclass Model_pred:\n    'Class for prediction with multiple models'\n    def __init__(self, models, use_tta=True, batch_size=32):\n        self.models = models\n        self.bs = batch_size\n        self.tfms = [tta.HorizontalFlip()] if use_tta else [] #, tta.VerticalFlip()]  \n        \n    def predict(self, ds):\n        #rasterio cannot be used with multiple workers\n        dl = DataLoader(ds, self.bs, num_workers=0, shuffle=False, pin_memory=True)\n        \n        # Create zero arrays\n        pred = np.zeros(ds.out_data_shape, dtype='float32')\n        merge_map = np.zeros(ds.out_data_shape, dtype='float32')\n        \n        # Gaussian weights\n        gw_numpy = _get_gaussian(ds.output_shape)\n        gw = torch.from_numpy(gw_numpy).to(device)\n        \n        with torch.no_grad():\n            for images, idxs in tqdm(iter(dl), total=len(dl)):\n                if ((idxs>=0).sum() > 0): #exclude empty images\n                    images = images[idxs>=0].to(device)\n                    idxs = idxs[idxs>=0]\n                    merger = tta.Merger()\n                    for t in tta.Compose(self.tfms):\n                        aug_images = t.augment_image(images)\n                        model_merger = tta.Merger()\n                        for model in tqdm(self.models):\n                            out = model(aug_images)\n                            out = F.softmax(out, dim=1)\n                            model_merger.append(out)\n                            torch.cuda.empty_cache()\n                        out = t.deaugment_mask(model_merger.result())\n                        merger.append(out)\n            \n                    # Apply gaussian weigthing\n                    batch_smx = merger.result()*gw.view(1,1,*gw.shape)\n                    batch_smx = [x for x in batch_smx.permute(0,2,3,1).cpu().numpy()]\n\n                    for smx, idx in zip(batch_smx, idxs):\n                        slcs = ds.out_slices[idx]\n                        # Only using positive class here\n                        pred[slcs] += smx[...,1]\n                        merge_map[slcs] += gw_numpy\n\n        pred /= merge_map\n        return pred","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_path = Path('../input/hubmap-single-fold-models-b1b5')\n    \n    # zoom factor (e.g., 3 means downscaling from 1536 to 512)\n    scale = 3 \n    # tile shift for prediction\n    shift = 0.8 \n    tile_shape = (512, 512)\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 = 4\n    \n    # test time augmentation\n    tta = True\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":"def 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":"\nmodels_path_list = [\n    '../input/hubmap-efficientnet-b4/model-fold-4.h5',\n    '../input/hubmap-efficientnet-b4/model-fold-0.h5',\n    '../input/hubmap-efficientnet-b4/model-fold-3.h5',\n    '../input/hubmap-efficientnet-b6-pseudo/model-fold-2.h5',\n    '../input/hubmap-efficientnet-b6-pseudo/model-fold-1.h5',\n    #'../input/hubmap-efficientnetb7-pseudo-labelled/model-fold-2.h5',\n]\n\nfold_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))\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"deepflash_models=[]\n\n'''\nfor fold in [1,2]:\n        \n        ###################### efficient net b4 #########################\n        model = smp.Unet(encoder_name='efficientnet-b4', \n                         encoder_weights=cfg.encoder_weights, \n                         in_channels=cfg.in_channels, \n                         classes=cfg.classes)\n        model, stats = load_model_weights(model, '../input/hubmap-deepflash-efficientnetb4/'+f'unet_efficientnet-b4_{fold}.pth')\n        \n        \n        if torch.cuda.is_available():  model.cuda()\n        \n        deepflash_models.append(model)\n        del model\n        gc.collect()\n        \n        model = smp.Unet(encoder_name='efficientnet-b3', \n                         encoder_weights=cfg.encoder_weights, \n                         in_channels=cfg.in_channels, \n                         classes=cfg.classes)\n        model, stats = load_model_weights(model, '../input/hubmap-deepflash-efficientnet-b3/'+f'unet_efficientnet-b3_{fold}.pth')\n        \n        if torch.cuda.is_available():  model.cuda()\n        \n        deepflash_models.append(model)\n        del model\n        gc.collect()\n'''      \n#         res_id=fold+2\n#         model = smp.Unet(encoder_name='resnet50', \n#                          encoder_weights=cfg.encoder_weights, \n#                          in_channels=cfg.in_channels, \n#                          classes=cfg.classes)\n#         model, stats = load_model_weights(model, '../input/hubmap-effiecient-sampling-deepflash-resnet/'+f'unet_resnet50_{res_id}.pth')\n#         if torch.cuda.is_available():  model.cuda()\n        \n#         deepflash_models.append(model)\n#         del model\n#         gc.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')\n\nif torch.cuda.is_available(): model.cuda()\n\ndeepflash_models.append(model)\ndel model\ngc.collect()\n\n# model = smp.Unet(encoder_name='efficientnet-b7', \n#                  encoder_weights=cfg.encoder_weights, \n#                  in_channels=cfg.in_channels, \n#                  classes=cfg.classes)\n# model, stats = load_model_weights(model, '../input/hubmap-deepflash-effiicentnetb7/'+f'unet_efficientnet-b7.pth')\n\n# if torch.cuda.is_available(): model.cuda()\n\n# deepflash_models.append(model)\n# del model\n# gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Sample submissions for ids\nimport 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[1:2, :]\nelse:\n    df_sample=df_sample\n\nprint(len(deepflash_models))\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Prediction","metadata":{}},{"cell_type":"code","source":"identity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n\nnames,predicts = [],[]\nfor idx,row in tqdm(df_sample.iterrows(),total=len(df_sample)):\n    \n    print(f'###### File {idx} ######')\n    f = cfg.data_path/'test'/f'{idx}.tiff'\n    ds = HubmapDataset(f, stats, scale=cfg.scale, shift=cfg.shift, output_shape=cfg.tile_shape)\n    \n    print('Predicting...')   \n    \n    dl = DataLoader(ds, 1, num_workers=0, shuffle=False, pin_memory=True)\n    \n    preds = np.zeros(ds.out_data_shape, dtype=np.uint8)\n    print(preds.shape)\n    with torch.no_grad():\n        for images, idxs in tqdm(iter(dl), total=len(dl)):\n            if ((idxs>=0).sum() > 0): #exclude empty images\n                images = images[idxs>=0].to(device)\n                idxs = idxs[idxs>=0]\n                merger = tta.Merger()\n                for t in tta.Compose([tta.HorizontalFlip(), tta.VerticalFlip()]):\n                    aug_images = t.augment_image(images)\n                    model_merger = tta.Merger()\n                    for model in deepflash_models:\n                        out = model(aug_images)\n                        #print(out.shape)\n                        out = F.softmax(out, dim=1)\n                        model_merger.append(out)\n                        torch.cuda.empty_cache()\n                    out = t.deaugment_mask(model_merger.result()).detach().cpu().numpy()/len(deepflash_models)\n                    preds[ds.out_slices[idxs]]=(out[0, 1, :, :] > 0.3).astype(np.uint8)\n                   \n    del ds, dl, images, merger, aug_images, model, out, model_merger\n    gc.collect()\n    torch.cuda.empty_cache()\n    \n    dataset_tf = rasterio.open(f, transform = identity) \n    slices_tf = make_grid_tf(dataset_tf.shape, window=WINDOW, min_overlap=MIN_OVERLAP)\n    print(dataset_tf.shape)\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    preds_tf = np.zeros(dataset_tf.shape, dtype=np.uint8)\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)\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                    K.clear_session()\n\n                pred_temp = pred_temp/len(fold_models)\n\n\n            pred_temp = cv2.resize((pred_temp).astype('uint8'), (WINDOW, WINDOW))\n            preds_tf[x1:x2,y1:y2] += pred_temp\n            \n            del image, hsv, s_th, p_th, pred_temp\n            gc.collect()\n\n    \n    del EMPTY, slices_tf, dataset_tf, fold_model\n    gc.collect()\n    K.clear_session()\n    \n    preds = cv2.resize(preds (preds_tf.shape[0], preds_tf.shape[1]))\n    preds = preds + preds_tf\n    preds = (preds >= VOTERS).astype(np.uint8)\n    del preds_tf\n    gc.collect()\n    \n#     plt.hist(preds)\n#     plt.show()\n\n    rle = rle_encode_less_memory(preds)\n    names.append(idx)\n    predicts.append(rle)\n    \n    del preds, rle\n    gc.collect()\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'id':names,'predicted':predicts})\ndf.to_csv('submission.csv',index=False)\ndf.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}