{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\nimport os\n\nimport glob\nfrom tqdm.notebook import tqdm\nimport shutil\nfrom pathlib import Path\nfrom openslide import OpenSlide\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage.measure import block_reduce\nimport cv2.cv2 as cv2\nfrom sys import getsizeof\nimport collections\nimport pickle","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nphysical_devices = tf.config.list_physical_devices('GPU')\ntf.config.experimental.set_memory_growth(physical_devices[0], True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = '../input/prostate-cancer-grade-assessment'\ndf_train = pd.read_csv(os.path.join(data_dir, 'train.csv'))\ndf_test = pd.read_csv(os.path.join(data_dir, 'test.csv'))\ndf_sub = pd.read_csv(os.path.join(data_dir, 'sample_submission.csv'))\n\nmodel_dir = '../input/panda-public-models'\nimage_folder = os.path.join(data_dir, 'test_images')\nis_test = os.path.exists(image_folder)  # IF test_images is not exists, we will use some train images.\nimage_folder = image_folder if is_test else os.path.join(data_dir, 'train_images')\n\ndf = df_test if is_test else df_train.loc[:10]\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def roll(a,      # ND array\n         window,      # rolling 2D window array\n         dx=1,   # horizontal step, abscissa, number of columns\n         dy=1):  # vertical step, ordinate, number of rows\n    shape = (\n        ((a.shape[0] - window[0]) // dy + 1),\n        ((a.shape[1] - window[1]) // dx + 1),\n        *window,\n    )\n    strides = (\n        a.strides[0] * dy,\n        a.strides[1] * dx,\n        a.strides[0],\n        a.strides[1],\n    )\n    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)\n\n\ndef reduce(image, block, stride):\n    return roll(image, (block, block), stride, stride).mean((-1,-2))\n\ndef tile_filter(tile, provider):\n    if provider[:3] == 'rad':\n        if ((1.0 - tile / 255).mean(2) < 0.1).astype(np.float).mean() > 0.4:\n            return False\n        if ((tile / 255).mean(2) < 0.4).astype(np.float).mean() > 0.2:\n            return False\n    else:\n        if ((1.0 - tile / 255).mean(2) < 0.2).astype(np.float).mean() > 0.3:\n            return False\n    return True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.models import load_model\nfrom tensorflow.keras.applications.resnet50 import preprocess_input\nfrom tensorflow.keras.preprocessing.image import img_to_array, load_img\nmodel = load_model('../input/tilenet/checkout_panda_v3.hdf5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"output = [model.output, model.layers[-2].output]\nK = output[1].shape[1:]\nfrom tensorflow.keras.models import Model\nmmodel = Model(inputs=model.input, outputs = output)\nprint(f'# features: {K}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## WSI Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"from typing import Type\nfrom torch import nn, Tensor\n\n\nclass Conv2dBlock(nn.Module):\n    def __init__(self, in_channels: int, out_channels: int, layers: int = 2,\n                 activation: Type[nn.Module] = nn.ReLU, kernel_size: int = 3,\n                 batchnorm: bool = True):\n        super(Conv2dBlock, self).__init__()\n        self.blocks = []\n\n        for i in range(layers):\n            block_layers = [\n                nn.Conv2d((in_channels if i == 0 else out_channels),\n                          out_channels, kernel_size, padding=kernel_size // 2),\n                activation()\n            ]\n            if batchnorm:\n                block_layers.insert(1, nn.BatchNorm2d(out_channels))\n            block = nn.Sequential(*block_layers)\n            self.blocks.append(block)\n\n        self.main = nn.Sequential(*self.blocks)\n\n    def forward(self, inputs: Tensor) -> Tensor:\n        \"\"\"Run layer\"\"\"\n        return self.main(inputs)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from torch.nn.functional import pad\nimport torch\nfrom typing import Union, Tuple, Iterable\n\nclass PandaModel(nn.Module):\n    def __init__(self, n_channels: int, n_classes: Union[int, Tuple[int, ...]]):\n        super().__init__()\n        self.encoder = nn.Sequential(\n            Conv2dBlock(n_channels, 1024),\n         #   nn.MaxPool2d(2),\n            Conv2dBlock(1024, 512),\n         #   nn.MaxPool2d(2),\n            Conv2dBlock(512, 256),\n         #   nn.MaxPool2d(2),\n            Conv2dBlock(256, 128),\n        )\n\n        def make_head(n: int):\n            return nn.Sequential(\n                nn.Linear(128, 32),\n                nn.ReLU(),\n                nn.Linear(32, 16),\n                nn.ReLU(),\n                nn.Linear(16, n)\n            )\n\n        if isinstance(n_classes, Iterable):\n            self.heads = nn.ModuleList([make_head(n) for n in n_classes])\n        else:\n            self.heads = nn.ModuleList([make_head(n_classes)])\n\n    def forward(self, x):\n        d = self.encoder(x).detach()\n        d = d.max(3).values.max(2).values.detach()\n        res = [h(d).detach() for h in self.heads]\n        return res if len(res) > 1 else res[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"wsi_model = PandaModel(2048, 6)\nwsi_model.load_state_dict(torch.load('../input/wsinet/best_qp2_723.pth', map_location='cpu')['model'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Normalization"},{"metadata":{"trusted":true},"cell_type":"code","source":"def color_deconvolution(im_rgb, w, I_0=None):\n    if np.linalg.norm(w[:, 2]) <= 1e-16:\n        wc = complement_stain_matrix(w)\n    else:\n        wc = w\n\n    wc = normalize(wc)\n    Q = np.linalg.inv(wc)\n    m = convert_image_to_matrix(im_rgb)[:3]\n    sda_fwd = rgb_to_sda(m, I_0)\n    sda_deconv = np.dot(Q, sda_fwd)\n    sda_inv = sda_to_rgb(sda_deconv, 255 if I_0 is not None else None)\n\n    # reshape output\n    StainsFloat = convert_matrix_to_image(sda_inv, im_rgb.shape)\n\n    # transform type\n    Stains = StainsFloat.clip(0, 255).astype(np.uint8)\n\n    # return\n    Unmixed = collections.namedtuple('Unmixed',\n                                     ['Stains', 'StainsFloat', 'Wc'])\n    Output = Unmixed(Stains, StainsFloat, wc)\n\n    return Output","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stain_color_map = {\n    'hematoxylin': [0.65, 0.70, 0.29],\n    'eosin':       [0.07, 0.99, 0.11],\n    'dab':         [0.27, 0.57, 0.78],\n    'null':        [0.0, 0.0, 0.0]\n}\n\ndef stain_unmixing_routine(\n        im_rgb, stains=None, stain_unmixing_method='macenko_pca',\n        stain_unmixing_params=None, mask_out=None):\n    stains = ['hematoxylin', 'eosin'] if stains is None else stains\n    stain_unmixing_params = {} if stain_unmixing_params is None else stain_unmixing_params\n\n    stain_unmixing_method = stain_unmixing_method.lower()\n\n    if stain_unmixing_method == 'macenko_pca':\n        stain_deconvolution = rgb_separate_stains_macenko_pca\n        stain_unmixing_params['I_0'] = None\n        stain_unmixing_params['mask_out'] = mask_out\n\n    elif stain_unmixing_method == 'xu_snmf':\n        stain_deconvolution = rgb_separate_stains_xu_snmf\n        stain_unmixing_params['I_0'] = None\n        assert mask_out is None, \"Masking is not yet implemented in xu_snmf.\"\n\n    else:\n        raise ValueError(\"Unknown/Unimplemented deconvolution method.\")\n\n    # get W_source\n    W_source = stain_deconvolution(im_rgb, **stain_unmixing_params)\n\n    # If Macenco method, reorder channels in W_target and W_source as desired.\n    # This is actually a necessary step in macenko's method since we're\n    # not guaranteed the order of the different stains.\n    if stain_unmixing_method == 'macenko_pca':\n        W_source = _reorder_stains(W_source, stains=stains)\n\n    return W_source\n\n\ndef color_deconvolution_routine(\n        im_rgb, W_source=None, mask_out=None, **kwargs):\n    # get W_source if not provided\n    if W_source is None:\n        W_source = stain_unmixing_routine(im_rgb, mask_out=mask_out, **kwargs)\n\n    # deconvolve\n    Stains, StainsFloat, wc = color_deconvolution(im_rgb, w=W_source, I_0=None)\n\n    # mask out (keep in mind, image is inverted)\n    if mask_out is not None:\n        for i in range(3):\n            Stains[..., i][mask_out] = 255\n            StainsFloat[..., i][mask_out] = 255.\n\n    return Stains, StainsFloat, wc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def find_stain_index(reference, w):\n    dot_products = np.dot(reference, w)\n    return np.argmax(np.abs(dot_products))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _reorder_stains(W, stains=None):\n    stains = ['hematoxylin', 'eosin'] if stains is None else stains\n\n    assert len(stains) == 2, \"Only two-stain matrices are supported for now.\"\n\n    def _get_channel_order(W):\n        first = find_stain_index(stain_color_map[stains[0]], W)\n        second = 1 - first\n        # If 2 stains, third \"stain\" is cross product of 1st 2 channels\n        # calculated using complement_stain_matrix()\n        third = 2\n        return first, second, third\n\n    def _ordered_stack(mat, order):\n        return np.stack([mat[..., j] for j in order], -1)\n\n    return _ordered_stack(W, _get_channel_order(W))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def convert_image_to_matrix(im):\n    if im.ndim == 2:\n        return im\n\n    return im.reshape((-1, im.shape[-1])).T\n\ndef convert_matrix_to_image(m, shape):\n    if len(shape) == 2:\n        return m\n\n    return m.T.reshape(shape[:-1] + (m.shape[0],))\n\nimport numpy\n\ndef get_principal_components(m):\n    return numpy.linalg.svd(m.astype(float), full_matrices=False)[0].astype(m.dtype)\n\n\ndef magnitude(m):\n    return numpy.sqrt((m ** 2).sum(0))\n\n\ndef normalize(m):\n    return m / magnitude(m)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def rgb_to_sda(im_rgb, I_0, allow_negatives=False):\n    is_matrix = im_rgb.ndim == 2\n    if is_matrix:\n        im_rgb = im_rgb.T\n\n    if I_0 is None:  # rgb_to_od compatibility\n        im_rgb = im_rgb.astype(float) + 1\n        I_0 = 256\n\n    if not allow_negatives:\n        im_rgb = np.minimum(im_rgb, I_0)\n\n    im_sda = -np.log(im_rgb/(1.*I_0)) * 255/np.log(I_0)\n    return im_sda.T if is_matrix else im_sda\n\ndef sda_to_rgb(im_sda, I_0):\n    is_matrix = im_sda.ndim == 2\n    if is_matrix:\n        im_sda = im_sda.T\n\n    od = I_0 is None\n    if od:  # od_to_rgb compatibility\n        I_0 = 256\n\n    im_rgb = I_0 ** (1 - im_sda / 255.)\n    return (im_rgb.T if is_matrix else im_rgb) - od\n\ndef exclude_nonfinite(m):\n    return m[:, np.isfinite(m).all(axis=0)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def color_convolution(im_stains, w, I_0=None):\n    m = convert_image_to_matrix(im_stains)\n    sda_fwd = rgb_to_sda(m, 255 if I_0 is not None else None,allow_negatives=True)\n    sda_conv = np.dot(w, sda_fwd)\n    sda_inv = sda_to_rgb(sda_conv, I_0)\n\n    # reshape output, transform type\n    im_rgb = (convert_matrix_to_image(sda_inv, im_stains.shape).clip(0, 255).astype(np.uint8))\n\n    return im_rgb","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def complement_stain_matrix(w):\n    stain0 = w[:, 0]\n    stain1 = w[:, 1]\n    stain2 = np.cross(stain0, stain1)\n    # Normalize new vector to have unit norm\n    return np.array([stain0, stain1, stain2 / np.linalg.norm(stain2)]).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def rgb_separate_stains_macenko_pca(im_rgb, I_0, *args, **kwargs):\n    im_sda = rgb_to_sda(im_rgb, I_0)\n    return separate_stains_macenko_pca(im_sda, *args, **kwargs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def separate_stains_macenko_pca(\n        im_sda, minimum_magnitude=16, min_angle_percentile=0.01,\n        max_angle_percentile=0.99, mask_out=None):\n    m = convert_image_to_matrix(im_sda)\n\n    # mask out irrelevant values\n    if mask_out is not None:\n        keep_mask = numpy.equal(mask_out[..., None], False)\n        keep_mask = numpy.tile(keep_mask, (1, 1, 3))\n        keep_mask = convert_image_to_matrix(keep_mask)\n        m = m[:, keep_mask.all(axis=0)]\n\n    # get rid of NANs and infinities\n    m = exclude_nonfinite(m)\n\n    # Principal components matrix\n    pcs = get_principal_components(m)\n    # Input pixels projected into the PCA plane\n    proj = pcs.T[:-1].dot(m)\n    # Pixels above the magnitude threshold\n    filt = proj[:, magnitude(proj) > minimum_magnitude]\n    # The \"angles\"\n    angles = _get_angles(filt)\n\n    # The stain vectors\n\n    def get_percentile_vector(p):\n        return pcs[:, :-1].dot(filt[:, argpercentile(angles, p)])\n\n    min_v = get_percentile_vector(min_angle_percentile)\n    max_v = get_percentile_vector(max_angle_percentile)\n\n    # The stain matrix\n    w = complement_stain_matrix(normalize(numpy.array([min_v, max_v]).T))\n    return w\n\n\ndef _get_angles(m):\n    m = normalize(m)\n    return (1 - m[1]) * numpy.sign(m[0])\n\ndef argpercentile(arr, p):\n    \"\"\"Calculate index in arr of element nearest the pth percentile.\"\"\"\n    # Index corresponding to percentile\n    i = int(p * arr.size + 0.5)\n    return numpy.argpartition(arr, i)[i]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stain_unmixing_routine_params = {\n    'stains': ['hematoxylin', 'eosin'],\n    'stain_unmixing_method': 'macenko_pca',\n}\n\ndef dbn(im_src, W_source=None, W_target=None, im_target=None,stains=None, mask_out=None, stain_unmixing_routine_params=None):\n    stains = ['hematoxylin', 'eosin'] if stains is None else stains\n    stain_unmixing_routine_params = (\n        {} if stain_unmixing_routine_params is None else\n        stain_unmixing_routine_params)\n    for k in ['W_source', 'mask_out']:\n        assert k not in stain_unmixing_routine_params.keys(), \\\n            \"%s must be provided as a separate parameter.\" % k\n\n    # find stains matrix from source image\n    stain_unmixing_routine_params['stains'] = stains\n    _, StainsFloat, _ = color_deconvolution_routine(\n        im_src, W_source=W_source, mask_out=mask_out,\n        **stain_unmixing_routine_params)\n\n    # Get W_target\n\n    if all(j is None for j in [W_target, im_target]):\n        # Normalize to 'ideal' stain matrix if none is provided\n        W_target = np.array(\n            [stain_color_map[stains[0]], stain_color_map[stains[1]]]).T\n        W_target = complement_stain_matrix(W_target)\n\n    elif im_target is not None:\n        # Get W_target from target image\n        W_target = stain_unmixing_routine(\n            im_target, **stain_unmixing_routine_params)\n\n    # Convolve source image StainsFloat with W_target\n    im_src_normalized = color_convolution(StainsFloat, W_target)\n\n    # return masked values using unnormalized image\n    if mask_out is not None:\n        keep_mask = np.not_equal(mask_out, True)\n        for i in range(3):\n            original = im_src[:, :, i].copy()\n            new = im_src_normalized[:, :, i].copy()\n            original[keep_mask] = 0\n            new[mask_out] = 0\n            im_src_normalized[:, :, i] = new + original\n\n    return im_src_normalized","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Predict"},{"metadata":{"trusted":true},"cell_type":"code","source":"LEVEL = 1\nTILE_SIZE = 256\n\nres_isup = np.zeros(len(df), dtype=np.int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"wsi_model = wsi_model.cuda()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import skimage.io\nbi = None\ns = None\ntiles = None\nfor i, row in tqdm(df.iterrows(), total=len(df)):\n    try:\n        image_id = row['image_id']\n        slide = skimage.io.MultiImage(f'{image_folder}/{image_id}.tiff', conserve_memory=False)[1]#OpenSlide(f'{image_folder}/{image_id}.tiff')\n        size = slide.shape[:2]\n        #size = slide.level_dimensions[2]\n        #downscale = slide.level_downsamples[2] / slide.level_downsamples[LEVEL]\n        #pos_downscale = slide.level_downsamples[LEVEL]\n\n        ys = np.arange(0,size[0],TILE_SIZE)\n        xs = np.arange(0,size[1],TILE_SIZE)\n        ty = len(ys)\n        tx = len(xs)\n        ys, xs = np.meshgrid(ys,xs)\n        ys = ys.reshape(-1)\n        xs = xs.reshape(-1)\n        feat_image = np.zeros((ty, tx, K[0]))\n        \n        tss = []\n        for y, x in zip(ys, xs):\n            tile = slide[y:y+TILE_SIZE, x:x+TILE_SIZE]\n            if not tile_filter(tile, row['data_provider']):\n                continue\n\n            mask_out = ((1.0 - tile / 255).mean(2) <  0.2)\n            tissue_rgb_normalized_kar = dbn(tile,\n                                            stain_unmixing_routine_params=stain_unmixing_routine_params,\n                                            mask_out=mask_out)\n            p, f = mmodel(preprocess_input(tissue_rgb_normalized_kar[None]))\n            feat_image[y // TILE_SIZE, x // TILE_SIZE] = f[0]\n            tss.append(tile)\n            \n        if bi is None or (bi == 0).mean() > (feat_image == 0).mean():\n            bi = feat_image\n            sss = slide\n            tiles = tss\n            print('CLICK')\n\n        with torch.no_grad():  \n            feats = torch.tensor(feat_image.astype(np.float32)).permute(2, 0, 1).cuda()\n            s = feats.shape[1:]\n            s = [16 - s[0], 16 - s[1]]\n            if s[0] < 0:\n                s[0] = 0\n            if s[1] < 0:\n                s[1] = 0\n            feats = pad(feats, [s[1] // 2, (s[1] + 1) // 2, s[0] // 2, (s[0] + 1) // 2, ]).detach()\n            res = wsi_model(feats[None]).cpu()\n            res = res.detach()[0].softmax(0)\n            print(res)\n            isup_score = res.argmax().item()\n            res_isup[i] = isup_score\n    except:\n        pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(sss)\ncv2.imwrite('3_wsi.png', cv2.cvtColor(sss, cv2.COLOR_RGB2BGR))\n\nk = 0\nfor t in tiles:\n    cv2.imwrite(f'3_tiles_{k}.png', cv2.cvtColor(t, cv2.COLOR_RGB2BGR))\n    k+=1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"m_ind = (bi == 0).mean((0,1))\nidxs = np.argsort(m_ind)[:3]\n\nk = 0\nfor i in idxs:\n    fig = plt.figure(frameon=False)\n    ax = plt.Axes(fig, [0., 0., 1., 1.])\n    ax.set_axis_off()\n    fig.add_axes(ax)\n    ax.imshow(bi[...,i], cmap='inferno')\n    fig.savefig(f'3_features_{k}.png')\n    k+=1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"df['isup_grade'] = np.array(res_isup, dtype=np.int)\ndf[['image_id', 'isup_grade']].to_csv('submission.csv', index=False)\nprint(df.head())\nprint()\nprint(df.isup_grade.value_counts())","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}