{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import cv2\nfrom tqdm import tqdm_notebook as tqdm\nimport fastai\nfrom fastai.vision import *\nimport os\nfrom mish_activation import *\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport skimage.io\nimport numpy as np\nimport pandas as pd\nsys.path.insert(0, '../input/semisupervised-imagenet-models/semi-supervised-ImageNet1K-models-master/')\nfrom hubconf import *\nimport math\nimport warnings","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#DATA = '../input/prostate-cancer-grade-assessment/train_images/'\n#TEST = '../input/prostate-cancer-grade-assessment/train.csv'\nDATA = '../input/prostate-cancer-grade-assessment/test_images'\nTEST = '../input/prostate-cancer-grade-assessment/test.csv'\nSAMPLE = '../input/prostate-cancer-grade-assessment/sample_submission.csv'\nMODELS = [f'../input/panda-init-class-model1/RNXT50_128krnew1_3featureB_{i}.pth' for i in range(4)] + \\\n         [f'../input/panda-init-class-model1/RNXT50_128krnew1_2featureB_{i}.pth' for i in range(4)] + \\\n         [f'../input/panda-init-class-model1/RNXT50_128krnew1_1featureB_{i}.pth' for i in range(4)]\n\nMODELS1 = [f'../input/panda-init-class-model1/RNXT50_128krnew_2_{i}.pth' for i in range(4)]\n\nsz = 128\nbs = 1\nN_max = 3*128\nN_max1 = 256\nnworkers = 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Modelm1(nn.Module):\n    def __init__(self, arch='resnext50_32x4d', n=6, pre=True):\n        super().__init__()\n        #m = torch.hub.load('facebookresearch/semi-supervised-ImageNet1K-models', arch)\n        m = _resnext(semi_supervised_model_urls[arch], Bottleneck, [3, 4, 6, 3], False, progress=False,\\\n                     groups=32,width_per_group=4)\n        self.enc = nn.Sequential(*list(m.children())[:-2])       \n        nc = list(m.children())[-1].in_features\n        #self.head = nn.Sequential(AdaptiveConcatPool2d(),Flatten(),nn.Linear(2*nc,512),Mish(),nn.BatchNorm1d(512),\n        #                          nn.Dropout(0.5),nn.Linear(512,n))\n        self.head = nn.Sequential(AdaptiveConcatPool2d(),Flatten(),\n                                  nn.Linear(2*nc,512),Mish(),nn.GroupNorm(32,512),\n                                  nn.Dropout(0.5),nn.Linear(512,n))\n        #self.head_s = nn.Conv2d(nc,7,1)\n        \n    def forward(self, x):\n        shape = x.shape\n        n = shape[1]\n        x = x.view(-1,shape[2],shape[3],shape[4])\n        x = self.enc(x)\n        shape = x.shape\n        x = x.view(-1,n,shape[1],shape[2],shape[3]).permute(0,2,1,3,4).contiguous()\\\n          .view(-1,shape[1],shape[2]*n,shape[3])\n        x = self.head(x)\n        return x[:,:1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AdaptiveConcatPool2dm1(Module):\n    \"Layer that concats `AdaptiveAvgPool2d` and `AdaptiveMaxPool2d`.\"\n    def __init__(self, sz:Optional[int]=None):\n        \"Output will be 2*sz or 2 if sz is None\"\n        self.output_size = sz or 1\n        self.ap = nn.AdaptiveAvgPool2d(self.output_size)\n        self.mp = nn.AdaptiveMaxPool2d(self.output_size)\n\n    def forward(self, x): return torch.cat([self.mp(x), self.ap(x)], 1)\n    #def forward(self, x): return torch.cat([self.ap(x),self.mp(x)], 1)#v1\n\ndef _resnext(url, block, layers, pretrained, progress, **kwargs):\n    model = ResNet(block, layers, **kwargs)\n    #state_dict = load_state_dict_from_url(url, progress=progress)\n    #model.load_state_dict(state_dict)\n    return model\n\nclass Model(nn.Module):\n    def __init__(self, arch='resnext50_32x4d', n=6, pre=True):\n        super().__init__()\n        #m = torch.hub.load('facebookresearch/semi-supervised-ImageNet1K-models', arch)\n        m = _resnext(semi_supervised_model_urls[arch], Bottleneck, [3, 4, 6, 3], False, progress=False,\\\n                     groups=32,width_per_group=4)\n        self.enc = nn.Sequential(*list(m.children())[:-2])       \n        nc = list(m.children())[-1].in_features\n        #self.head = nn.Sequential(AdaptiveConcatPool2d(),Flatten(),nn.Linear(2*nc,512),Mish(),nn.BatchNorm1d(512),\n        #                          nn.AdaptiveConcatPool2dm1(0.5),nn.Linear(512,n))\n        self.head = nn.Sequential(AdaptiveConcatPool2d(),Flatten(),\n                                  nn.Linear(2*nc,512),Mish(),nn.GroupNorm(32,512),\n                                  nn.Dropout(0.5),nn.Linear(512,n))\n        #self.head_s = nn.Conv2d(nc,7,1)\n        \n    def forward(self, x):\n        shape = x.shape\n        n = shape[1]\n        x = x.view(-1,shape[2],shape[3],shape[4])\n        x = self.enc(x)\n        shape = x.shape\n        x = x.view(-1,n,shape[1],shape[2],shape[3]).permute(0,2,1,3,4).contiguous()\\\n          .view(-1,shape[1],shape[2]*n,shape[3])\n        x = self.head(x)\n        return x[:,:1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AdaptiveConcatPool2dm(Module):\n    \"Layer that concats `AdaptiveAvgPool2d` and `AdaptiveMaxPool2d`.\"\n    def __init__(self, sz:Optional[int]=None):\n        \"Output will be 2*sz or 2 if sz is None\"\n        self.output_size = sz or 1\n        self.ap = nn.AdaptiveAvgPool2d(self.output_size)\n        self.mp = nn.AdaptiveMaxPool2d(self.output_size)\n\n    #def forward(self, x): return torch.cat([self.mp(x), self.ap(x)], 1)\n    def forward(self, x): return torch.cat([self.ap(x),self.mp(x)], 1)#v1\n\ndef _resnext(url, block, layers, pretrained, progress, **kwargs):\n    model = ResNet(block, layers, **kwargs)\n    #state_dict = load_state_dict_from_url(url, progress=progress)\n    #model.load_state_dict(state_dict)\n    return model\n\nclass Modelm(nn.Module):\n    def __init__(self, arch='resnext50_32x4d', n=6, pre=True):\n        super().__init__()\n        #m = torch.hub.load('facebookresearch/semi-supervised-ImageNet1K-models', arch)\n        m = _resnext(semi_supervised_model_urls[arch], Bottleneck, [3, 4, 6, 3], False, progress=False,\\\n                     groups=32,width_per_group=4)\n        self.enc = nn.Sequential(*list(m.children())[:-2])       \n        nc = list(m.children())[-1].in_features\n        #self.head = nn.Sequential(AdaptiveConcatPool2d(),Flatten(),nn.Linear(2*nc,512),Mish(),nn.BatchNorm1d(512),\n        #                          nn.Dropout(0.5),nn.Linear(512,n))\n        self.head = nn.Sequential(AdaptiveConcatPool2dm(),Flatten(),\n                                  nn.Linear(2*nc,512),Mish(),nn.GroupNorm(32,512),\n                                  nn.Dropout(0.5),nn.Linear(512,n))\n        #self.head_s = nn.Conv2d(nc,7,1)\n        \n    def forward(self, x):\n        shape = x.shape\n        n = shape[1]\n        x = x.view(-1,shape[2],shape[3],shape[4])\n        x = self.enc(x)\n        shape = x.shape\n        x = x.view(-1,n,shape[1],shape[2],shape[3]).permute(0,2,1,3,4).contiguous()\\\n          .view(-1,shape[1],shape[2]*n,shape[3])\n        x = self.head(x)\n        return x[:,:1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"models = []\nfor path in MODELS[:-4]:\n    state_dict = torch.load(path,map_location=torch.device('cpu'))\n    model = Model(n=1+10)\n    #model = Model()\n    model.load_state_dict(state_dict)\n    model.float()\n    model.eval()\n    model.cuda()\n    models.append(model)\nfor path in MODELS[-4:]:\n    state_dict = torch.load(path,map_location=torch.device('cpu'))\n    model = Modelm(n=1+10)\n    #model = Model()\n    model.load_state_dict(state_dict)\n    model.float()\n    model.eval()\n    model.cuda()\n    models.append(model)\n    \nmodels2 = []\nfor path in MODELS1:\n    state_dict = torch.load(path,map_location=torch.device('cpu'))\n    model = Modelm1(n=1+10)\n    #model = Model()\n    model.load_state_dict(state_dict)\n    model.float()\n    model.eval()\n    model.cuda()\n    models2.append(model)\n\ndel state_dict","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _mask_tissue(image, kernel_size=(7, 7), gray_threshold=220):\n    \"\"\"Masks tissue in image. Uses gray-scaled image, as well as\n    dilation kernels and 'gap filling'\n    \"\"\"\n    # Define elliptic kernel\n    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, kernel_size)\n    # Convert rgb to gray scale for easier masking\n    gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)\n    # Now mask the gray-scaled image (capturing tissue in biopsy)\n    mask = np.where(gray < gray_threshold, 1, 0).astype(np.uint8)\n    # Use dilation and findContours to fill in gaps/holes in masked tissue\n    mask = cv2.dilate(mask, kernel, iterations=1)\n    contour, _ = cv2.findContours(mask, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)\n    for cnt in contour:\n        cv2.drawContours(mask, [cnt], 0, 1, -1)\n    return mask\n\ndef _pad_image(image, pad_len, pad_val):\n    \"\"\"Pads inputted image, accepts both \n    2-d (mask) and 3-d (rgb image) arrays\n    \"\"\"\n    if image is None:\n        return None\n    elif image.ndim == 2:\n        return np.pad(\n            image, ((pad_len, pad_len), (pad_len, pad_len)), pad_val)\n    elif image.ndim == 3:\n        return np.pad(\n            image, ((pad_len, pad_len), (pad_len, pad_len), (0, 0)), pad_val)\n    return None\n\ndef _transpose_image(image):\n    \"\"\"Inputs an image and transposes it, accepts \n    both 2-d (mask) and 3-d (rgb image) arrays\n    \"\"\"\n    if image is None:\n        return None\n    elif image.ndim == 2:\n        return np.transpose(image, (1, 0)).copy()\n    elif image.ndim == 3:\n        return np.transpose(image, (1, 0, 2)).copy()\n    return None\n\ndef _get_tissue_parts_indices(tissue, min_consec_info):\n    \"\"\"If there are multiple tissue parts in 'tissue', 'tissue' will be \n    split it. Each tissue part will be taken care of separately, and if \n    the tissue part is less than min_consec_info, it's considered to small\n    and won't be returned.\n    \"\"\"\n    split_points = np.where(np.diff(tissue) != 1)[0]+1\n    tissue_parts = np.split(tissue, split_points)\n    return [\n        tp for tp in tissue_parts if len(tp) >= min_consec_info\n    ]\n\ndef _get_tissue_subparts_coords(subtissue, patch_size, min_decimal_keep):\n    \"\"\"Inputs a tissue part resulting from '_get_tissue_parts_indices'.\n    This tissue part is divided into N subparts and returned.\n    Argument min_decimal_keep basically decides if we should squeeze in the \n    N subparts in an area bigger than the sum of the N subparts or not.\n    \"\"\"\n    start, end = subtissue[0], subtissue[-1]\n    num_subparts = (end-start)/patch_size\n    if num_subparts % 1 < min_decimal_keep and num_subparts >= 1:\n        num_subparts = math.floor(num_subparts)\n    else:\n        num_subparts = math.ceil(num_subparts)\n\n    excess = (num_subparts*patch_size) - (end-start)\n    shift = excess // 2\n\n    return [\n        i * patch_size + start - shift \n        for i in range(num_subparts)\n    ]\n\ndef _eval_and_append_xy_coords(coords,\n                               image, \n                               mask, \n                               patch_size, \n                               x, y, \n                               min_patch_info,\n                               transposed,\n                               precompute):\n    \"\"\"Based on computed x and y coordinates of patch: \n    slices out patch from original image, flattens it,\n    preprocesses it, and finally evaluates its mask.\n    If patch contains more info than min_patch_info,\n    the patch coordinates are kept, along with a \n    value 'val1' that estimates how much information\n    there is in the patch.\n    \"\"\"\n    patch_1d = (\n        image[y: y+patch_size, x:x+patch_size, :]\n        .mean(axis=2)\n        .reshape(-1)\n    )\n    idx_tissue = np.where(patch_1d <= 210)[0]\n    idx_black = np.where(patch_1d < 5)[0]\n    idx_background = np.where(patch_1d > 210)[0]\n\n    if len(idx_tissue) > 0:\n        patch_1d[idx_black] = 210\n        patch_1d[idx_background] = 210\n        val1 = int(patch_1d.mean())\n        val2 = mask[y:y+patch_size, x:x+patch_size].mean()\n        if val2 > min_patch_info:\n            if precompute:\n                if transposed:\n                    coords = np.concatenate([\n                        coords, [[val1, x-patch_size, y-patch_size]]\n                    ])\n                else:\n                    coords = np.concatenate([\n                        coords, [[val1, y-patch_size, x-patch_size]]\n                    ])\n            else:\n                coords = np.concatenate([\n                    coords, [[val1, y, x]]\n                ])\n               \n    return coords\n\ndef compute_coords(image,\n                   patch_size=256,\n                   precompute=False,\n                   min_patch_info=0.35,\n                   min_axis_info=0.35,\n                   min_consec_axis_info=0.35,\n                   min_decimal_keep=0.7):\n\n    \"\"\"\n    Input:\n        image : 3-d np.ndarray\n        patch_size : size of patches/tiles, will be of \n            size (patch_size x patch_size x 3)\n        precompute : If True, only coordinates will be returned,\n            these coordinates match the inputted 'original' image.\n            If False, both an image and coordinates will be returned,\n            the coordinates does not match the inputted image but the\n            image that it is returned with.\n        min_patch_info : Minimum required information in patch\n            (see '_eval_and_append_xy_coords')\n        min_axis_info : Minimum fraction of on-bits in x/y dimension to be \n            considered enough information. For x, this would be fraction of \n            on-bits in x-dimension of a y:y+patch_size slice. For y, this would \n            be the fraction of on-bits for the whole image in y-dimension\n        min_consec_axis_info : Minimum consecutive x/y on-bits from 'min_axis_info'\n            (see '_get_tissue_parts_indices')\n        min_decimal_keep : Threshold for decimal point for removing \"excessive\" patch\n            (see '_get_tissue_subparts_coords')\n    \n    Output:\n        image [only if precompute is True] : similar to input image, but fits \n            to the computed coordinates\n        coords : the coordinates that will be used to compute the patches later on\n    \"\"\"\n    \n    \n    if type(image) != np.ndarray:\n        # if image is a Tensor\n        image = image.numpy()\n    \n    # masked tissue will be used to compute the coordinates\n    mask = _mask_tissue(image)\n\n    # initialize coordinate accumulator\n    coords = np.zeros([0, 3], dtype=int)\n\n    # pad image and mask to make sure no tissue is potentially missed out\n    image = _pad_image(image, patch_size, 'maximum')\n    mask = _pad_image(mask, patch_size, 'minimum')\n    \n    y_sum = mask.sum(axis=1)\n    x_sum = mask.sum(axis=0)\n    # if on bits in x_sum is greater than in y_sum, the tissue is\n    # likely aligned horizontally. The algorithm works better if\n    # the image is aligned vertically, thus the image will be transposed\n    if len(np.where(x_sum > 0)[0]) > len(np.where(y_sum > 0)[0]):\n        image = _transpose_image(image)\n        mask = _transpose_image(mask)\n        y_sum, _ = x_sum, y_sum\n        transposed = True\n    else:\n        transposed = False\n    \n    # where y_sum is more than the minimum number of on-bits\n    y_tissue = np.where(y_sum >= (patch_size*min_axis_info))[0]\n    \n    if len(y_tissue) < 1:\n        warnings.warn(\"Not enough tissue in image (y-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    y_tissue_parts_indices = _get_tissue_parts_indices(\n        y_tissue, patch_size*min_consec_axis_info)\n    \n    if len(y_tissue_parts_indices) < 1: \n        warnings.warn(\"Not enough tissue in image (y-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    # loop over the tissues in y-dimension\n    for yidx in y_tissue_parts_indices:\n        y_tissue_subparts_coords = _get_tissue_subparts_coords(\n            yidx, patch_size, min_decimal_keep)\n        \n        for y in y_tissue_subparts_coords:\n            # in y_slice, where x_slice_sum is more than the minimum number of on-bits\n            x_slice_sum = mask[y:y+patch_size, :].sum(axis=0)\n            x_tissue = np.where(x_slice_sum >= (patch_size*min_axis_info))[0]\n            \n            x_tissue_parts_indices = _get_tissue_parts_indices(\n                x_tissue, patch_size*min_consec_axis_info)\n            \n            # loop over tissues in x-dimension (inside y_slice 'y:y+patch_size')\n            for xidx in x_tissue_parts_indices:\n                x_tissue_subparts_coords = _get_tissue_subparts_coords(\n                    xidx, patch_size, min_decimal_keep)\n                \n                for x in x_tissue_subparts_coords:\n                    coords = _eval_and_append_xy_coords(\n                        coords, image, mask, patch_size, x, y, \n                        min_patch_info, transposed, precompute\n                    )     \n    \n    if len(coords) < 1:\n        warnings.warn(\"Not enough tissue in image (x-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    if precompute: return coords\n    else: return image, coords\n    \ndef tile(img):\n    coords = compute_coords(img,\n                            patch_size=sz,\n                            precompute=True,\n                            min_patch_info=0.35,\n                            min_axis_info=0.35,\n                            min_consec_axis_info=0.35,\n                            min_decimal_keep=0.7)\n    tiles = []\n    for (v, y, x) in coords:\n        p = img[y:y+sz,x:x+sz]\n        if p.shape[0] != sz or p.shape[1] != sz: continue\n        tiles.append(p)\n    if len(tiles) == 0: return None\n    tiles = [x for y, x in sorted(zip([t.sum() for t in tiles], tiles), key=lambda pair: pair[0])]\n    if len(tiles) > N_max: tiles = tiles[:N_max]\n    t1 = np.stack(tiles)\n    if N_max1 > len(tiles): \n        for i in range(N_max1 - len(tiles)): tiles.append(np.full((sz,sz,3), 255, dtype=np.uint8))\n    if len(tiles) > N_max1: tiles = tiles[:N_max1]\n    t2 = np.stack(tiles)\n    return t1,t2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean = torch.Tensor([1.0-0.85506157, 1.0-0.7035249, 1.0-0.80203127])\nstd = torch.Tensor([0.40011922, 0.52504386, 0.42675745])\n\nclass PandaDataset(Dataset):\n    def __init__(self, path, test):\n        self.path = path\n        self.names = list(pd.read_csv(test).image_id)\n\n    def __len__(self):\n        return len(self.names)\n\n    def __getitem__(self, idx):\n        name = self.names[idx]\n        img = skimage.io.MultiImage(os.path.join(DATA,name+'.tiff'))[1]\n        t1,t2 = tile(img)\n        tiles = torch.Tensor((255 - t1)/255.0)\n        tiles = (tiles - mean)/std\n        tiles2 = torch.Tensor((255 - t2)/255.0)\n        tiles2 = (tiles2 - mean)/std\n        return tiles.permute(0,3,1,2),tiles2.permute(0,3,1,2), name","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df = pd.read_csv(SAMPLE)\nif os.path.exists(DATA):\n    ds = PandaDataset(DATA,TEST)\n    dl = DataLoader(ds, batch_size=bs, num_workers=nworkers, shuffle=False)\n    names,preds = [],[]\n\n    with torch.no_grad():\n        for x,x2,y in tqdm(dl):\n            x = x.cuda()\n            #dihedral TTA\n            x = torch.stack([x,x.flip(-1),x.flip(-2),x.flip(-1,-2),x.transpose(-1,-2),\\\n              x.transpose(-1,-2).flip(-1), x.transpose(-1,-2).flip(-2),\\\n              x.transpose(-1,-2).flip(-1,-2)],1)\n            x = x.view(8,-1,3,sz,sz)\n            p = [model(x[a].unsqueeze(0)) for model in models for a in range(x.shape[0])]\n            x2 = x2.cuda()\n            x2 = torch.stack([x2,x2.flip(-1),x2.flip(-2),x2.flip(-1,-2),x2.transpose(-1,-2),\\\n              x2.transpose(-1,-2).flip(-1), x2.transpose(-1,-2).flip(-2),\\\n              x2.transpose(-1,-2).flip(-1,-2)],1)\n            x2 = x2.view(8,-1,3,sz,sz)\n            p += [model(x2[a].unsqueeze(0)) for model in models2 for a in range(x2.shape[0])]\n            \n            p = torch.stack(p,1)\n            #p = p.view(bs,8*len(models),-1).mean(1).argmax(-1).cpu()\n            #print(6.0*torch.sigmoid(p.view(bs,8*len(models))))\n            p = p.view(bs,8*(len(models)+len(models2))).mean(1).cpu()\n            #print(6.0*torch.sigmoid(p))\n            p = torch.clamp((6.0*torch.sigmoid(p)).long().view(-1),0,5)\n            \n            names.append(y)\n            preds.append(p)\n    \n    names = np.concatenate(names)\n    preds = torch.cat(preds).numpy()\n    sub_df = pd.DataFrame({'image_id': names, 'isup_grade': preds})\n    sub_df.to_csv('submission.csv', index=False)\n    sub_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df.to_csv(\"submission.csv\", index=False)\nsub_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}