{"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":"# Submission notebook of HuBMAP - Pytorch smp Unet++ Inference\n\n\nTrain part is here.\n\nhttps://www.kaggle.com/nayuts/hubmap-pytorch-smp-unet","metadata":{}},{"cell_type":"markdown","source":"--------------------------","metadata":{}},{"cell_type":"markdown","source":"I create train and inference notebook for Unet++ of [segmentation_models.pytorch](https://github.com/qubvel/segmentation_models.pytorch).\n\nI use only pytorch for framework.\n\nThe accuracy of the model is going to be tuned and improved in the future.\n\nI published this notebook for our reference as an example implementation using only pytorch.","metadata":{"papermill":{"duration":0.018407,"end_time":"2020-11-29T03:22:16.658942","exception":false,"start_time":"2020-11-29T03:22:16.640535","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"I refered following two great notebooks for training and inference.\n\n- https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub\n\n- https://www.kaggle.com/curiosity806/hubmap-use-catalyst-smp-and-albumentations\n\nAnd refered following great notebook for dice loss.\n\n- https://www.kaggle.com/bigironsphere/loss-function-library-keras-pytorch\n\nKindly upvote and appreciate the original work.","metadata":{}},{"cell_type":"markdown","source":"## Load libraries","metadata":{"papermill":{"duration":0.016657,"end_time":"2020-11-29T03:22:16.726414","exception":false,"start_time":"2020-11-29T03:22:16.709757","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import gc\nimport os\nimport random\nimport sys\nimport time\nimport warnings\nwarnings.simplefilter(\"ignore\")\n\n#import pdb\n#import zipfile\n#import pydicom\nfrom albumentations import *\nfrom albumentations.pytorch import ToTensor\nimport cv2\nfrom matplotlib import pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image, ImageFilter\nimport rasterio\nfrom rasterio.windows import Window\nfrom sklearn.model_selection import KFold\nimport tifffile as tiff\nimport torch\nimport torch.backends.cudnn as cudnn\nimport torch.nn as nn\nfrom torch.nn import functional as F\nimport torch.optim as optim\nfrom torch.optim.lr_scheduler import ReduceLROnPlateau\nfrom torch.utils.data import DataLoader, Dataset, sampler\nfrom tqdm import tqdm_notebook as tqdm","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":4.529346,"end_time":"2020-11-29T03:22:41.66512","exception":false,"start_time":"2020-11-29T03:22:37.135774","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To use segmentation_models_pytorch in offline environment, I loaded modules I cloned and uploaded to dataset.","metadata":{}},{"cell_type":"code","source":"#https://www.kaggle.com/hfutybx/unet-densenet121-lung-of-segmentation/data\n\nsys.path.append('../input/efficientnetpytorchaug252020/EfficientNet-PyTorch-master')\nsys.path.append('../input/pretrainedmodels/pretrainedmodels-0.7.4/')\nsys.path.append('../input/pytorchimagemodelsoct302020/pytorch-image-models-master')\nsys.path.append('../input/segmentation-models-pytorch0-1-2/segmentation_models.pytorch-master')\nimport segmentation_models_pytorch as smp","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##  Set parameters","metadata":{"papermill":{"duration":0.027134,"end_time":"2020-11-29T03:22:41.720618","exception":false,"start_time":"2020-11-29T03:22:41.693484","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def set_seed(seed=46):\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    np.random.seed(seed)\n    random.seed(seed)\n    torch.backends.cudnn.deterministic = True\nset_seed(46)","metadata":{"papermill":{"duration":0.038393,"end_time":"2020-11-29T03:22:41.786754","exception":false,"start_time":"2020-11-29T03:22:41.748361","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To save time, the number of \"folds\" is smaller. \nThese days, depending on the model, I think it's more common to use 5 ~ 10.","metadata":{"papermill":{"duration":0.027878,"end_time":"2020-11-29T03:22:41.843067","exception":false,"start_time":"2020-11-29T03:22:41.815189","status":"completed"},"tags":[]}},{"cell_type":"code","source":"fold = 0\nnfolds = 5\nreduce = 4\nsz = 256\n\nBATCH_SIZE = 16\nDEVICE = ('cuda' if torch.cuda.is_available() else 'cpu')\nNUM_WORKERS = 4\nSEED = 2020\nTH = 0.50  #threshold for positive predictions\n\nDATA = '../input/hubmap-kidney-segmentation/test/'\nLABELS = '../input/hubmap-kidney-segmentation/train.csv'\ndf_sample = pd.read_csv('../input/hubmap-kidney-segmentation/sample_submission.csv')","metadata":{"papermill":{"duration":0.395969,"end_time":"2020-11-29T03:22:42.268205","exception":false,"start_time":"2020-11-29T03:22:41.872236","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Util functions","metadata":{"papermill":{"duration":0.029423,"end_time":"2020-11-29T03:22:42.329459","exception":false,"start_time":"2020-11-29T03:22:42.300036","status":"completed"},"tags":[]}},{"cell_type":"code","source":"#https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\n#with bug fix\ndef rle_encode_less_memory(img):\n    #watch out for the bug\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":{"papermill":{"duration":0.038315,"end_time":"2020-11-29T03:22:42.395843","exception":false,"start_time":"2020-11-29T03:22:42.357528","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset","metadata":{"papermill":{"duration":0.028236,"end_time":"2020-11-29T03:22:42.451856","exception":false,"start_time":"2020-11-29T03:22:42.42362","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def get_aug(p=1.0):\n    return Compose([\n        HorizontalFlip(),\n        VerticalFlip(),\n        RandomRotate90(),\n        ShiftScaleRotate(shift_limit=0.0625, scale_limit=0.2, rotate_limit=15, p=0.9, \n                         border_mode=cv2.BORDER_REFLECT),\n        OneOf([\n            OpticalDistortion(p=0.3),\n            GridDistortion(p=.1),\n            IAAPiecewiseAffine(p=0.3),\n        ], p=0.3),\n        OneOf([\n            HueSaturationValue(10,15,10),\n            CLAHE(clip_limit=2),\n            RandomBrightnessContrast(),            \n        ], p=0.3),\n    ], p=p)","metadata":{"papermill":{"duration":0.039197,"end_time":"2020-11-29T03:22:42.596336","exception":false,"start_time":"2020-11-29T03:22:42.557139","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model","metadata":{"papermill":{"duration":0.105237,"end_time":"2020-11-29T03:22:57.078897","exception":false,"start_time":"2020-11-29T03:22:56.97366","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def get_UnetPlusPlus():\n    model =  smp.UnetPlusPlus(\n                 encoder_name='efficientnet-b3',\n                 encoder_weights=None,\n                 in_channels=3,\n                 classes=1)\n    return model","metadata":{"papermill":{"duration":0.110928,"end_time":"2020-11-29T03:22:57.293516","exception":false,"start_time":"2020-11-29T03:22:57.182588","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference\n\nBefore notebook version8, I created simple dataset class which gets list of tiled images, and return. But new test data, there are big image data, the memory allocates. For this reason, I am reusing the modules from [HuBMAP Pytorch/fast.ai starter sub](https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub).","metadata":{"papermill":{"duration":0.074016,"end_time":"2020-11-29T03:54:36.68492","exception":false,"start_time":"2020-11-29T03:54:36.610904","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# https://www.kaggle.com/iafoss/256x256-images\nmean = np.array([0.63701495, 0.4709702,  0.6817423] )\nstd = np.array([0.15978882, 0.2245109, 0.14173926])\n\ndef img2tensor(img,dtype:np.dtype=np.float32):\n    if img.ndim==2 : img = np.expand_dims(img,2)\n    img = np.transpose(img,(2,0,1))\n    return torch.from_numpy(img.astype(dtype, copy=False))\n\nclass HuBMAPTestDataset(Dataset):\n    def __init__(self, idx, sz=sz, reduce=reduce):\n        self.data = rasterio.open(os.path.join(DATA,idx+'.tiff'), transform = identity,\n                                 num_threads='all_cpus')\n        # some images have issues with their format \n        # and must be saved correctly before reading with rasterio\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        self.shape = self.data.shape\n        self.reduce = reduce\n        self.sz = reduce*sz\n        self.pad0 = (self.sz - self.shape[0]%self.sz)%self.sz\n        self.pad1 = (self.sz - self.shape[1]%self.sz)%self.sz\n        self.n0max = (self.shape[0] + self.pad0)//self.sz\n        self.n1max = (self.shape[1] + self.pad1)//self.sz\n        \n    def __len__(self):\n        return self.n0max*self.n1max\n    \n    def __getitem__(self, idx):\n        # the code below may be a little bit difficult to understand,\n        # but the thing it does is mapping the original image to\n        # tiles created with adding padding, as done in\n        # https://www.kaggle.com/iafoss/256x256-images ,\n        # and then the tiles are loaded with rasterio\n        # n0,n1 - are the x and y index of the tile (idx = n0*self.n1max + n1)\n        n0,n1 = idx//self.n1max, idx%self.n1max\n        # x0,y0 - are the coordinates of the lower left corner of the tile in the image\n        # negative numbers correspond to padding (which must not be loaded)\n        x0,y0 = -self.pad0//2 + n0*self.sz, -self.pad1//2 + n1*self.sz\n        # make sure that the region to read is within the image\n        p00,p01 = max(0,x0), min(x0+self.sz,self.shape[0])\n        p10,p11 = max(0,y0), min(y0+self.sz,self.shape[1])\n        img = np.zeros((self.sz,self.sz,3),np.uint8)\n        # mapping the loade region to the tile\n        if self.data.count == 3:\n            img[(p00-x0):(p01-x0),(p10-y0):(p11-y0)] = np.moveaxis(self.data.read([1,2,3],\n                window=Window.from_slices((p00,p01),(p10,p11))), 0, -1)\n        else:\n            for i,layer in enumerate(self.layers):\n                img[(p00-x0):(p01-x0),(p10-y0):(p11-y0),i] =\\\n                  layer.read(1,window=Window.from_slices((p00,p01),(p10,p11)))\n        \n        if self.reduce != 1:\n            img = cv2.resize(img,(self.sz//reduce,self.sz//reduce),\n                             interpolation = cv2.INTER_AREA)\n        #check for empty imges\n        hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n        h,s,v = cv2.split(hsv)\n        if (s>s_th).sum() <= p_th or img.sum() <= p_th:\n            #images with -1 will be skipped\n            return img2tensor((img/255.0 - mean)/std), -1\n        else: return img2tensor((img/255.0 - mean)/std), idx","metadata":{"papermill":{"duration":0.085068,"end_time":"2020-11-29T03:54:36.84431","exception":false,"start_time":"2020-11-29T03:54:36.759242","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!ls ../input/hubmap-pytorch-smp-unet/","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models = []\n\nfor fold in range(nfolds):\n\n    model = get_UnetPlusPlus().to(DEVICE)\n    model.load_state_dict(torch.load(f\"../input/hubmap-pytorch-smp-unet/FOLD{fold}_.pth\"))\n    models.append(model)","metadata":{"papermill":{"duration":1.250809,"end_time":"2020-11-29T03:54:38.169407","exception":false,"start_time":"2020-11-29T03:54:36.918598","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#iterator like wrapper that returns predicted masks\nclass Model_pred:\n    def __init__(self, models, dl, tta:bool=True, half:bool=False):\n        self.models = models\n        self.dl = dl\n        self.tta = tta\n        self.half = half\n        \n    def __iter__(self):\n        count=0\n        with torch.no_grad():\n            for x,y in iter(self.dl):\n                if ((y>=0).sum() > 0): #exclude empty images\n                    x = x[y>=0].to(DEVICE)\n                    y = y[y>=0]\n                    if self.half: x = x.half()\n                    py = None\n                    for model in self.models:\n                        p = model(x)\n                        p = torch.sigmoid(p).detach()\n                        if py is None: py = p\n                        else: py += p\n                    if self.tta:\n                        #x,y,xy flips as TTA\n                        flips = [[-1],[-2],[-2,-1]]\n                        for f in flips:\n                            xf = torch.flip(x,f)\n                            for model in self.models:\n                                p = model(xf)\n                                p = torch.flip(p,f)\n                                py += torch.sigmoid(p).detach()\n                        py /= (1+len(flips))        \n                    py /= len(self.models)\n\n                    py = F.upsample(py, scale_factor=reduce, mode=\"bilinear\")\n                    py = py.permute(0,2,3,1).float().cpu()\n                    \n                    batch_size = len(py)\n                    for i in range(batch_size):\n                        yield py[i],y[i]\n                        count += 1\n                    \n    def __len__(self):\n        return len(self.dl.dataset)","metadata":{"papermill":{"duration":0.094455,"end_time":"2020-11-29T03:54:38.339388","exception":false,"start_time":"2020-11-29T03:54:38.244933","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s_th = 40  #saturation blancking threshold\np_th = 1000*(sz//256)**2 #threshold for the minimum number of pixels\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n\nnames,preds = [],[]\nprint(df_sample)\nfor idx,row in tqdm(df_sample.iterrows(),total=len(df_sample)):\n    idx = row['id']\n    ds = HuBMAPTestDataset(idx)\n    #rasterio cannot be used with multiple workers\n    dl = DataLoader(ds,BATCH_SIZE,num_workers=0,shuffle=False,pin_memory=True)\n    mp = Model_pred(models,dl)\n    #generate masks\n    mask = torch.zeros(len(ds),ds.sz,ds.sz,dtype=torch.int8)\n    for p,i in iter(mp): mask[i.item()] = p.squeeze(-1) > TH\n    \n    #reshape tiled masks into a single mask and crop padding\n    mask = mask.view(ds.n0max,ds.n1max,ds.sz,ds.sz).\\\n        permute(0,2,1,3).reshape(ds.n0max*ds.sz,ds.n1max*ds.sz)\n    mask = mask[ds.pad0//2:-(ds.pad0-ds.pad0//2) if ds.pad0 > 0 else ds.n0max*ds.sz,\n        ds.pad1//2:-(ds.pad1-ds.pad1//2) if ds.pad1 > 0 else ds.n1max*ds.sz]\n    \n    #convert to rle\n    #https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\n    rle = rle_encode_less_memory(mask.numpy())\n    names.append(idx)\n    preds.append(rle)\n    del mask, ds, dl\n    gc.collect()","metadata":{"papermill":{"duration":590.629497,"end_time":"2020-11-29T04:04:29.322675","exception":false,"start_time":"2020-11-29T03:54:38.693178","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'id':names,'predicted':preds})\ndf.to_csv('submission.csv',index=False)","metadata":{"papermill":{"duration":1.363307,"end_time":"2020-11-29T04:04:30.761583","exception":false,"start_time":"2020-11-29T04:04:29.398276","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]}]}