{"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.4","metadata":{"id":"OcsetTMwKXqC"}},{"cell_type":"markdown","source":"### Installation and package loading","metadata":{}},{"cell_type":"code","source":"# Install deepflash2 and dependencies\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","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from geojson import Point, Feature, FeatureCollection, dump\nimport shapely.wkt\nimport shapely.geometry\nfrom shapely.geometry import MultiPolygon, Polygon\nimport matplotlib.pyplot as plt\nimport cv2, torch, tifffile, pandas as pd, gc\nfrom fastai.vision.all import *","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import rasterio","metadata":{"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\n# def read_img(self:BaseDataset, *args, **kwargs):\n#     image = tifffile.imread(args[0])\n#     if len(image.shape) == 5:\n#         image = image.squeeze().transpose(1, 2, 0)\n#     elif image.shape[0] == 3:\n#         image = image.transpose(1, 2, 0)\n#     return image\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\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/models/unet_efficientnet-b4_1.pth' #unet_efficientnet-b4_1\n    #../input/models/unet_efficientnet-b0_newlabel_eff0_1.pth\n    \n    # deepflash2 dataset (https://matjesg.github.io/deepflash2/data.html#TileDataset)\n    scale = 1.5 # zoom facor (zoom out) 3\n    tile_shape = (512, 512)\n    padding = (100,100) # Border overlap for prediction\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 = 16\n    \n    # prediction threshold\n    threshold = 0.5\n    \ncfg = CONFIG()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Sample submissions for ids\ndf_sample = pd.read_csv(cfg.data_path/'sample_submission.csv',  index_col='id')\n\n# Model (see https://github.com/qubvel/segmentation_models.pytorch)\nmodel = smp.Unet(encoder_name=cfg.encoder_name, \n                 encoder_weights=cfg.encoder_weights, \n                 in_channels=cfg.in_channels, \n                 classes=cfg.classes)\nmodel.segmentation_head[0] = nn.Conv2d(16, 1, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))\nmodel, stats = load_model_weights(model, cfg.model_file)\nbatch_tfms = [Normalize.from_stats(*stats)]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Prediction","metadata":{}},{"cell_type":"code","source":"import deepflash2.tta as tta","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cell\n@patch\ndef predict_tiles(self:Learner, ds_idx=1, dl=None, path=None, mc_dropout=False, n_times=1, use_tta=False,\n                       tta_merge='mean', tta_tfms=None, uncertainty_estimates=True, energy_T=1):\n    \"Make predictions and reconstruct tiles, optional with dropout and/or tta applied.\"\n\n    if dl is None: dl = self.dls[ds_idx].new(shuffled=False, drop_last=False)\n    assert isinstance(dl.dataset, TileDataset), \"Provide dataloader containing a TileDataset\"\n    if use_tta: tfms = tta_tfms or [tta.HorizontalFlip(), tta.Rotate90(angles=[90,180,270])]\n    else: tfms=[]\n\n    self.model.eval()\n    if mc_dropout: self.apply_dropout()\n\n    store = str(path) if path else zarr.storage.TempStore()\n    root = zarr.group(store=store, overwrite=True)\n    g_smx, g_seg, g_std, g_eng  = root.create_groups('smx', 'seg', 'std', 'energy')\n\n    i = 0\n    last_file = None\n    for data in progress_bar(dl, leave=False):\n        if isinstance(data, TensorImage): images = data\n        else: images, _, _ = data\n        m_smx = tta.Merger()\n        m_energy = tta.Merger()\n        out_list_smx = []\n        for t in tta.Compose(tfms):\n            for _ in range(n_times):\n                aug_images = t.augment_image(images)\n                with torch.no_grad():\n                    out = self.model(aug_images)\n                out = t.deaugment_mask(out)\n                if dl.padding[0]!= images.shape[-1]-out.shape[-1]:\n                    padding = ((images.shape[-1]-out.shape[-1]-dl.padding[0])//2,)*4\n                    out = F.pad(out, padding)\n                m_smx.append(sigmoid(out,))#change softmax(out,) to sigmoid(out,)\n                if uncertainty_estimates:\n                    e = (energy_T*torch.logsumexp(out/energy_T, dim=1)) #negative energy score\n                    m_energy.append(e)\n\n        ll = []\n        ll.append([x for x in m_smx.result().permute(0,2,3,1).cpu().numpy()])\n        if uncertainty_estimates:\n            ll.append([x for x in torch.mean(m_smx.result('std'), 1).cpu().numpy()])\n            ll.append([x for x in m_energy.result().cpu().numpy()])\n        for j, preds in enumerate(zip(*ll)):\n            if len(preds)==3: smx,std,eng = preds\n            else: smx = preds[0]\n            idx = i+j\n            f = dl.files[dl.image_indices[idx]]\n            outShape = dl.image_shapes[idx]\n            outSlice = dl.out_slices[idx]\n            inSlice = dl.in_slices[idx]\n            if last_file!=f:\n                z_smx = g_smx.zeros(f.name, shape=(*outShape, dl.c), dtype='float32')\n                z_seg = g_seg.zeros(f.name, shape=outShape, dtype='uint8')\n                z_std = g_std.zeros(f.name, shape=outShape, dtype='float32')\n                z_eng = g_eng.zeros(f.name, shape=outShape, dtype='float32')\n                last_file = f\n            z_smx[outSlice] = smx[inSlice]\n            z_seg[outSlice] = np.argmax(smx, axis=-1)[inSlice]\n            if uncertainty_estimates:\n                z_std[outSlice] = std[inSlice]\n                z_eng[outSlice] = eng[inSlice]\n        i += dl.bs\n\n    return g_smx, g_seg, g_std, g_eng\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"names,preds = [],[]\n\n\nfor idx, _ in df_sample.iterrows():\n    print(f'###### File {idx} ######')\n    f = cfg.data_path/'test'/f'{idx}.tiff'\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    # 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    if torch.cuda.is_available(): dls.cuda(), model.cuda()\n    learn = Learner(dls, model, loss_func='')\n    \n    # Predict tiles, see https://matjesg.github.io/deepflash2/learner.html#Learner.predict_tiles\n    print('Prediction')\n    dl_test = dls.train\n    dl_test.c = 1 \n    res = learn.predict_tiles(dl=dl_test, path='/kaggle/temp/', use_tta=False, uncertainty_estimates=False)\n    \n    th = cfg.threshold\n    \n    # Load mask from softmax prediction > threshold\n    msk = (res[0][f.name][:,:,0]>=th).astype(np.uint8) #[:,:,0] change shape (7996, 15907,1) to (7996, 15907)\n    print('Rezising')\n    msk = cv2.resize(msk, (shape[1], shape[0]))\n    rle = rle_encode_less_memory(msk)\n    names.append(idx)\n    preds.append(rle)\n    \n    # Plot Result\n    print('Plotting')\n    fig, ax = plt.subplots(figsize=(15,15))\n    ax.imshow(cv2.resize(msk, (1024, 1024)))\n    plt.show()\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')]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":[]},{"cell_type":"markdown","source":"# substitue my json","metadata":{}},{"cell_type":"code","source":"# def polygons_to_mask(polygons):\n#     img_mask = np.zeros(im_size, np.uint8)\n#     if not polygons:\n#         return img_mask\n#     int_coords = lambda x: np.array(x).round().astype(np.int64)#int32\n#     exteriors = [int_coords(poly.exterior.coords) for poly in polygons]\n#     interiors = [int_coords(pi.coords) for poly in polygons\n#                  for pi in poly.interiors]\n#     cv2.fillPoly(img_mask, exteriors, 1)\n#     cv2.fillPoly(img_mask, interiors, 0)\n#     return img_mask","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# list_test = [\"d488c759a\"]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df_info = pd.read_csv('../input/hubmap-kidney-segmentation/HuBMAP-20-dataset_information.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# b = df_info.height_pixels.tolist()\n# a = df_info.width_pixels.tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# names2 = df_info.image_file.tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# list_name = []\n# for name in names2:\n#     name = name[:9]\n#     list_name.append(name)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# shapes = []\n# for x,y in zip(b,a):\n#     c = (x,y)\n#     shapes.append(c)\n\n# dic_shape2 = dict(zip(list_name,shapes))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for idx in list_test:\n#     with open(f'../input/handmodify1/modified02.json') as f:\n#       shapes = json.load(f)\n#     all_pol = []\n#     print(idx)\n#     for shape in shapes:\n#         shape = shape['geometry']['coordinates']\n#         shape_numpy = np.array(shape,dtype=object)\n#         shape_numpy = shape_numpy.squeeze()\n\n#         pol = Polygon(shape_numpy)\n#         all_pol.append(pol)\n#     all_polygons = MultiPolygon(all_pol)\n#     if not all_polygons.is_valid:\n#         all_polygons = all_polygons.buffer(0)\n#     # Sometimes buffer() converts a simple Multipolygon to just a Polygon,\n#     # need to keep it a Multi throughout\n#     if all_polygons.type == 'Polygon':\n#         all_polygons = MultiPolygon([all_polygons])\n#     im_size = dic_shape2[idx]\n#     mk = polygons_to_mask(all_polygons)\n#     rle = rle_encode_less_memory(mk)\n\n#     print(mk.shape)\n    \n#     #show plots\n#     fig, ax = plt.subplots(figsize=(15,15))\n#     ax.imshow(cv2.resize(mk, (1024, 1024)))\n#     plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df_sample.loc[idx] = rle","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Submission","metadata":{}},{"cell_type":"code","source":"# df = pd.DataFrame({'id':names,'predicted':preds})\n# df.to_csv('submission.csv',index=False)\n# df.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df_sample.to_csv('submission.csv')\n# display(df_sample)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}