{"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":"# Description\nThis kernel provides a starter Pytorch code for inference that performs dividing the images into tiles([based on this kernel](https://www.kaggle.com/iafoss/256x256-images)), selection of tiles with tissue, evaluation of the predictions of multiple models with TTA, combining the tile masks back into image level masks, and conversion into RLE. The inference is performed based on models trained in the [fast.ai starter kernel](https://www.kaggle.com/iafoss/hubmap-fast-ai-starter), provided by me. I hope it will help you to get started with this competition.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nimport tifffile as tiff\nimport cv2\nimport os\nimport gc\nfrom tqdm.notebook import tqdm\n\nfrom fastai.vision.all import *\nfrom torch.utils.data import Dataset, DataLoader\nimport rasterio\nfrom rasterio.windows import Window\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bs = 8#32\nsz = 2*512   #the size of tiles\noverlap = 2*2*64\nreduce = 2#4 #reduce the original images by 4 times \n\nTH = 0.43  #threshold for positive predictions\nDATA = '../input/hubmap-kidney-segmentation/test/'\nMODELS = [f'../input/hubmap-models/model0_r2_{i}.pth' for i in range(4)]\ndf_sample = pd.read_csv('../input/hubmap-kidney-segmentation/sample_submission.csv')\nNUM_WORKERS = 2\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data","metadata":{}},{"cell_type":"code","source":"#functions to convert encoding to mask and mask to encoding\ndef enc2mask(encs, shape):\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for m,enc in enumerate(encs):\n        if isinstance(enc,np.float) and np.isnan(enc): continue\n        s = enc.split()\n        for i in range(len(s)//2):\n            start = int(s[2*i]) - 1\n            length = int(s[2*i+1])\n            img[start:start+length] = 1 + m\n    return img.reshape(shape).T\n\ndef mask2enc(mask, n=1):\n    pixels = mask.T.flatten()\n    encs = []\n    for i in range(1,n+1):\n        p = (pixels == i).astype(np.int8)\n        if p.sum() == 0: encs.append(np.nan)\n        else:\n            p = np.concatenate([[0], p, [0]])\n            runs = np.where(p[1:] != p[:-1])[0] + 1\n            runs[1::2] -= runs[::2]\n            encs.append(' '.join(str(x) for x in runs))\n    return encs\n\n#https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\ndef rle_encode_less_memory(img):\n    '''\n    img: numpy array, 1 - mask, 0 - background\n    Returns run length as string formated\n    This simplified method requires first and last pixel to be zero\n    '''\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":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/iafoss/256x256-images\nmean = np.array([0.65459856,0.48386562,0.69428385])\nstd = np.array([0.15167958,0.23584107,0.13146145])\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\n# https://www.kaggle.com/iafoss/256x256-images\ns_th = 10  #saturation blancking threshold\np_th = 1000*(sz//256)**2 #threshold for the minimum number of pixel\nnames,preds = [],[]\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n\nclass HuBMAPDataset(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 - overlap)\n        self.sz0 = reduce*sz\n        self.pad0 = (self.sz - (self.shape[0] - reduce*overlap)%self.sz)%self.sz\n        self.pad1 = (self.sz - (self.shape[1] - reduce*overlap)%self.sz)%self.sz\n        if self.pad0 < reduce*overlap: self.pad0 += self.sz\n        if self.pad1 < reduce*overlap: self.pad1 += self.sz\n        \n        self.n0 = (self.shape[0] - reduce*overlap + self.pad0)//self.sz\n        self.n1 = (self.shape[1] - reduce*overlap + self.pad1)//self.sz\n        \n    def __len__(self):\n        return self.n0*self.n1\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 dividing into tiles and mapping \n        # the image to tiles built with padding\n        n0,n1 = idx//self.n1, idx%self.n1\n        # x0,y0 - are the coordinates of the lower left corner of the tile\n        x0,y0 = -self.pad0//2 + n0*self.sz, -self.pad1//2 + n1*self.sz\n        # making sure that the region to read is within the image\n        p00,p01 = max(0,x0), min(x0+self.sz0,self.shape[0])\n        p10,p11 = max(0,y0), min(y0+self.sz0,self.shape[1])\n        img = np.zeros((self.sz0,self.sz0,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.sz0//reduce,self.sz0//reduce),interpolation = cv2.INTER_AREA)\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            return img2tensor((img/255.0 - mean)/std), -1\n        else: return img2tensor((img/255.0 - mean)/std), idx","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model","metadata":{}},{"cell_type":"code","source":"class FPN(nn.Module):\n    def __init__(self, input_channels:list, output_channels:list):\n        super().__init__()\n        self.convs = nn.ModuleList(\n            [nn.Sequential(nn.Conv2d(in_ch, out_ch*2, kernel_size=3, padding=1),\n             nn.ReLU(inplace=True), nn.BatchNorm2d(out_ch*2),\n             nn.Conv2d(out_ch*2, out_ch, kernel_size=3, padding=1))\n            for in_ch, out_ch in zip(input_channels, output_channels)])\n        \n    def forward(self, xs:list, last_layer):\n        hcs = [F.interpolate(c(x),scale_factor=2**(len(self.convs)-i),mode='bilinear') \n               for i,(c,x) in enumerate(zip(self.convs, xs))]\n        hcs.append(last_layer)\n        return torch.cat(hcs, dim=1)\n\nclass UnetBlock(Module):\n    def __init__(self, up_in_c:int, x_in_c:int, nf:int=None, blur:bool=False,\n                 self_attention:bool=False, **kwargs):\n        super().__init__()\n        self.shuf = PixelShuffle_ICNR(up_in_c, up_in_c//2, blur=blur, **kwargs)\n        self.bn = nn.BatchNorm2d(x_in_c)\n        ni = up_in_c//2 + x_in_c\n        nf = nf if nf is not None else max(up_in_c//2,32)\n        self.conv1 = ConvLayer(ni, nf, norm_type=None, **kwargs)\n        self.conv2 = ConvLayer(nf, nf, norm_type=None,\n            xtra=SelfAttention(nf) if self_attention else None, **kwargs)\n        self.relu = nn.ReLU(inplace=True)\n\n    def forward(self, up_in:Tensor, left_in:Tensor) -> Tensor:\n        s = left_in\n        up_out = self.shuf(up_in)\n        cat_x = self.relu(torch.cat([up_out, self.bn(s)], dim=1))\n        return self.conv2(self.conv1(cat_x))\n        \nclass _ASPPModule(nn.Module):\n    def __init__(self, inplanes, planes, kernel_size, padding, dilation, groups=1):\n        super().__init__()\n        self.atrous_conv = nn.Conv2d(inplanes, planes, kernel_size=kernel_size,\n                stride=1, padding=padding, dilation=dilation, bias=False, groups=groups)\n        self.bn = nn.BatchNorm2d(planes)\n        self.relu = nn.ReLU()\n\n        self._init_weight()\n\n    def forward(self, x):\n        x = self.atrous_conv(x)\n        x = self.bn(x)\n\n        return self.relu(x)\n\n    def _init_weight(self):\n        for m in self.modules():\n            if isinstance(m, nn.Conv2d):\n                torch.nn.init.kaiming_normal_(m.weight)\n            elif isinstance(m, nn.BatchNorm2d):\n                m.weight.data.fill_(1)\n                m.bias.data.zero_()\n\nclass ASPP(nn.Module):\n    def __init__(self, inplanes=512, mid_c=256, dilations=[6, 12, 18, 24], out_c=None):\n        super().__init__()\n        self.aspps = [_ASPPModule(inplanes, mid_c, 1, padding=0, dilation=1)] + \\\n            [_ASPPModule(inplanes, mid_c, 3, padding=d, dilation=d,groups=4) for d in dilations]\n        self.aspps = nn.ModuleList(self.aspps)\n        self.global_pool = nn.Sequential(nn.AdaptiveMaxPool2d((1, 1)),\n                        nn.Conv2d(inplanes, mid_c, 1, stride=1, bias=False),\n                        nn.BatchNorm2d(mid_c), nn.ReLU())\n        out_c = out_c if out_c is not None else mid_c\n        self.out_conv = nn.Sequential(nn.Conv2d(mid_c*(2+len(dilations)), out_c, 1, bias=False),\n                                    nn.BatchNorm2d(out_c), nn.ReLU(inplace=True))\n        self.conv1 = nn.Conv2d(mid_c*(2+len(dilations)), out_c, 1, bias=False)\n        self._init_weight()\n\n    def forward(self, x):\n        x0 = self.global_pool(x)\n        xs = [aspp(x) for aspp in self.aspps]\n        x0 = F.interpolate(x0, size=xs[0].size()[2:], mode='bilinear', align_corners=True)\n        x = torch.cat([x0] + xs, dim=1)\n        return self.out_conv(x)\n    \n    def _init_weight(self):\n        for m in self.modules():\n            if isinstance(m, nn.Conv2d):\n                torch.nn.init.kaiming_normal_(m.weight)\n            elif isinstance(m, nn.BatchNorm2d):\n                m.weight.data.fill_(1)\n                m.bias.data.zero_()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from torchvision.models.resnet import ResNet, Bottleneck\nclass UneXt50(nn.Module):\n    def __init__(self, stride=1, **kwargs):\n        super().__init__()\n        #encoder\n        m = ResNet(Bottleneck, [3, 4, 6, 3], groups=32, width_per_group=4)\n        #m = torch.hub.load('facebookresearch/semi-supervised-ImageNet1K-models',\n        #                   'resnext50_32x4d_ssl')\n        self.enc0 = nn.Sequential(m.conv1, m.bn1, nn.ReLU(inplace=True))\n        self.enc1 = nn.Sequential(nn.MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1),\n                            m.layer1) #256\n        self.enc2 = m.layer2 #512\n        self.enc3 = m.layer3 #1024\n        self.enc4 = m.layer4 #2048\n        #aspp with customized dilatations\n        self.aspp = ASPP(2048,256,out_c=512,dilations=[stride*1,stride*2,stride*3,stride*4])\n        self.drop_aspp = nn.Dropout2d(0.5)\n        #decoder\n        self.dec4 = UnetBlock(512,1024,256)\n        self.dec3 = UnetBlock(256,512,128)\n        self.dec2 = UnetBlock(128,256,64)\n        self.dec1 = UnetBlock(64,64,32)\n        self.fpn = FPN([512,256,128,64],[16]*4)\n        self.drop = nn.Dropout2d(0.1)\n        self.final_conv = ConvLayer(32+16*4, 1, ks=1, norm_type=None, act_cls=None)\n        \n    def forward(self, x):\n        enc0 = self.enc0(x)\n        enc1 = self.enc1(enc0)\n        enc2 = self.enc2(enc1)\n        enc3 = self.enc3(enc2)\n        enc4 = self.enc4(enc3)\n        enc5 = self.aspp(enc4)\n        dec3 = self.dec4(self.drop_aspp(enc5),enc3)\n        dec2 = self.dec3(dec3,enc2)\n        dec1 = self.dec2(dec2,enc1)\n        dec0 = self.dec1(dec1,enc0)\n        x = self.fpn([enc5, dec3, dec2, dec1], dec0)\n        x = self.final_conv(self.drop(x))\n        x = F.interpolate(x,scale_factor=2,mode='bilinear')\n        return x","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models = []\nfor path in MODELS:\n    state_dict = torch.load(path,map_location=torch.device('cpu'))\n    model = UneXt50(2)\n    model.load_state_dict(state_dict)\n    model.float()\n    model.eval()\n    model.to(device)\n    models.append(model)\n\ndel state_dict","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{}},{"cell_type":"code","source":"names,preds = [],[]\nfor idx,row in tqdm(df_sample.iterrows(),total=len(df_sample)):\n    idx = row['id']\n    ds = HuBMAPDataset(idx)\n    #rasterio cannot be used with multiple workers\n    dl = DataLoader(ds,bs,num_workers=0,shuffle=False,pin_memory=True)\n    mp = Model_pred(models,dl)\n    #generate masks\n    \n    mask = torch.zeros(ds.shape[0]+ds.pad0,ds.shape[1]+ds.pad1,dtype=torch.int8)\n    for p,id in iter(mp):\n        ##print(p.shape,mask.shape)\n        j,i = id%ds.n1,id//ds.n1\n        ix,iy = ds.sz*i, ds.sz*j\n        px0, px1 = reduce*overlap//2, -reduce*overlap//2\n        py0, py1 = reduce*overlap//2, -reduce*overlap//2\n        mask[(ix+px0):(ix+ds.sz0+px1),(iy+py0):(iy+ds.sz0+py1)] = \\\n          (p.squeeze(-1) > TH)[px0:(ds.sz0 + px1),py0:(ds.sz0 + py1)]\n\n    mask = mask[ds.pad0//2:-(ds.pad0-ds.pad0//2) if ds.pad0 > 0 else ds.shape[0],\n        ds.pad1//2:-(ds.pad1-ds.pad1//2) if ds.pad1 > 0 else ds.shape[1]]\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()\n    #break","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'id':names,'predicted':preds})\ndf.to_csv('submission0.csv',index=False)\ndf = df.set_index('id')\n#df.update(pd.read_csv('../input/d48-hand-labelled/d48_hand_labelled.csv').set_index('id'))\ndfr = pd.read_csv('../input/hubmap-models/submission_0938.csv')\ndfr = dfr.loc[dfr.id == 'd488c759a']\ndf.update(dfr.set_index('id'))\ndf = df.reset_index()\ndf.to_csv('submission.csv',index=False)","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0"},"execution_count":null,"outputs":[]}]}