{"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":"code","source":"DEBUG = False","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Initialization","metadata":{}},{"cell_type":"markdown","source":"### Packages","metadata":{}},{"cell_type":"code","source":"!mkdir -p /tmp/pip/cache/\n!cp ../input/segmentationmodelspytorch/segmentation_models/efficientnet_pytorch-0.6.3.xyz /tmp/pip/cache/efficientnet_pytorch-0.6.3.tar.gz\n!cp ../input/segmentationmodelspytorch/segmentation_models/pretrainedmodels-0.7.4.xyz /tmp/pip/cache/pretrainedmodels-0.7.4.tar.gz\n!cp ../input/segmentationmodelspytorch/segmentation_models/segmentation-models-pytorch-0.1.2.xyz /tmp/pip/cache/segmentation_models_pytorch-0.1.2.tar.gz\n!cp ../input/segmentationmodelspytorch/segmentation_models/timm-0.1.20-py3-none-any.whl /tmp/pip/cache/\n!cp ../input/segmentationmodelspytorch/segmentation_models/timm-0.2.1-py3-none-any.whl /tmp/pip/cache/\n!pip install --no-index --find-links /tmp/pip/cache/ efficientnet-pytorch\n!pip install --no-index --find-links /tmp/pip/cache/ segmentation-models-pytorch","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Imports","metadata":{}},{"cell_type":"code","source":"import os\nimport gc\nimport sys\nimport cv2\nimport json\nimport glob\nimport torch\nimport random\nimport rasterio\nimport numpy as np\nimport pandas as pd\n\nfrom tqdm.notebook import tqdm\nfrom collections import Counter\nfrom rasterio.windows import Window\nfrom matplotlib import pyplot as plt\nfrom torch.utils.data import DataLoader, Dataset","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Seeding","metadata":{}},{"cell_type":"code","source":"def seed_everything(seed):\n    \"\"\"\n    Seeds basic parameters for reproductibility of results.\n\n    Args:\n        seed (int): Number of the seed.\n    \"\"\"\n\n    random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = True\n    \n    \nSEED = 2021\nseed_everything(SEED)","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Params","metadata":{}},{"cell_type":"code","source":"CLASSES = [\"ftus\"]\n\nNUM_CLASSES = len(CLASSES)\n\nMEAN = np.array([0.66437738, 0.50478148, 0.70114894])\nSTD = np.array([0.15825711, 0.24371008, 0.13832686])\n\nDEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n\nIDENTITY = rasterio.Affine(1, 0, 0, 0, 1, 0)\nFLIPS = [[-1], [-2], [-2, -1]]\n\nDATA_PATH = '../input/hubmap-kidney-segmentation/'\nIMG_PATH = DATA_PATH + 'test/'","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data","metadata":{}},{"cell_type":"markdown","source":"### Load","metadata":{}},{"cell_type":"code","source":"import tifffile as tiff\n\n\ndef load_image(img_path, df_info, reduce_factor=1):\n    \"\"\"\n    Load image and make sure sizes matches df_info\n    \"\"\"\n    image_fname = img_path.rsplit(\"/\", -1)[-1]\n    \n    \n    W = int(df_info[df_info.image_file == image_fname][\"width_pixels\"])\n    H = int(df_info[df_info.image_file == image_fname][\"height_pixels\"])\n\n    img = tiff.imread(img_path).squeeze()\n\n    channel_pos = np.argwhere(np.array(img.shape) == 3)[0][0]\n    W_pos = np.argwhere(np.array(img.shape) == W)[0][0]\n    H_pos = np.argwhere(np.array(img.shape) == H)[0][0]\n\n    img = np.moveaxis(img, (H_pos, W_pos, channel_pos), (0, 1, 2))\n    \n    if reduce_factor > 1:\n        img = cv2.resize(\n            img,\n            (img.shape[1] // reduce_factor, img.shape[0] // reduce_factor),\n            interpolation=cv2.INTER_AREA,\n        )\n        \n    return img\n","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Transforms","metadata":{}},{"cell_type":"code","source":"import albumentations as albu\nfrom albumentations.pytorch import ToTensorV2\n\ndef HE_preprocess(augment=False, visualize=False, mean=MEAN, std=STD):\n    \"\"\"\n    Returns transformations for the H&E images.\n\n    Args:\n        augment (bool, optional): Whether to apply augmentations. Defaults to True.\n        visualize (bool, optional): Whether to use transforms for visualization. Defaults to False.\n        mean ([type], optional): Mean for normalization. Defaults to MEAN.\n        std ([type], optional): Standard deviation for normalization. Defaults to STD.\n\n    Returns:\n        albumentation transforms: transforms.\n    \"\"\"\n    if visualize:\n        normalizer = albu.Compose(\n            [albu.Normalize(mean=[0, 0, 0], std=[1, 1, 1]), ToTensorV2()], p=1\n        )\n    else:\n        normalizer = albu.Compose(\n            [albu.Normalize(mean=mean, std=std), ToTensorV2()], p=1\n        )\n    \n    if augment:\n        raise NotImplementedError\n\n    return normalizer\n","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### RLE encoding","metadata":{}},{"cell_type":"code","source":"def rle_encode_less_memory(img):\n    pixels = img.T.flatten()\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    return ' '.join(str(x) for x in runs)\n\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):\n            continue\n        enc_split = enc.split()\n        for i in range(len(enc_split) // 2):\n            start = int(enc_split[2 * i]) - 1\n            length = int(enc_split[2 * i + 1])\n            img[start: start + length] = 1 + m\n    return img.reshape(shape).T","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Dataset","metadata":{}},{"cell_type":"code","source":"class InferenceDataset(Dataset):\n    def __init__(\n        self,\n        original_img_path,\n        rle=None,\n        overlap_factor=1,\n        tile_size=256,\n        reduce_factor=4,\n        transforms=None,\n    ):\n        self.original_img = load_image(original_img_path, full_size=reduce_factor > 1)\n        self.orig_size = self.original_img.shape\n\n        # self.original_img = lab_normalization(self.original_img)\n\n        self.raw_tile_size = tile_size\n        self.reduce_factor = reduce_factor\n        self.tile_size = tile_size * reduce_factor\n\n        self.overlap_factor = overlap_factor\n\n        self.positions = self.get_positions()\n        self.transforms = transforms\n\n        if rle is not None:\n            self.mask = enc2mask(rle, (self.orig_size[1], self.orig_size[0])) > 0\n        else:\n            self.mask = None\n\n    def __len__(self):\n        return len(self.positions)\n\n    def get_positions(self):\n        top_x = np.arange(\n            0,\n            self.orig_size[0],  # +self.tile_size,\n            int(self.tile_size / self.overlap_factor),\n        )\n        top_y = np.arange(\n            0,\n            self.orig_size[1],  # +self.tile_size,\n            int(self.tile_size / self.overlap_factor),\n        )\n        starting_positions = []\n        for x in top_x:\n            right_space = self.orig_size[0] - (x + self.tile_size)\n            if right_space > 0:\n                boundaries_x = (x, x + self.tile_size)\n            else:\n                boundaries_x = (x + right_space, x + right_space + self.tile_size)\n\n            for y in top_y:\n                down_space = self.orig_size[1] - (y + self.tile_size)\n                if down_space > 0:\n                    boundaries_y = (y, y + self.tile_size)\n                else:\n                    boundaries_y = (y + down_space, y + down_space + self.tile_size)\n                starting_positions.append((boundaries_x, boundaries_y))\n\n        return starting_positions\n\n    def __getitem__(self, idx):\n        pos_x, pos_y = self.positions[idx]\n        img = self.original_img[pos_x[0]: pos_x[1], pos_y[0]: pos_y[1], :]\n\n        # img = lab_normalization(img)\n\n        # down scale to tile size\n        if self.reduce_factor > 1:\n            img = cv2.resize(\n                img, (self.raw_tile_size, self.raw_tile_size), interpolation=cv2.INTER_AREA\n            )\n\n        if self.transforms:\n            img = self.transforms(image=img)[\"image\"]\n\n        pos = np.array([pos_x[0], pos_x[1], pos_y[0], pos_y[1]])\n\n        return img, pos","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class InferenceEfficientDataset(InferenceDataset):\n    \"\"\"\n    Refs : \n    https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub/\n    https://www.kaggle.com/finlay/pytorch-fcn-resnet50-in-20-minute\n    \"\"\"\n    def __init__(\n        self,\n        original_img_path,\n        rle=None,\n        overlap_factor=1,\n        tile_size=256,\n        reduce_factor=4,\n        transforms=None,\n    ):\n            \n        self.raw_tile_size = tile_size\n        self.reduce_factor = reduce_factor\n        self.tile_size = tile_size * reduce_factor\n        \n        self.overlap_factor = overlap_factor\n        self.transforms = transforms\n\n        # Load image with rasterio        \n        self.original_img = rasterio.open(original_img_path, transform=IDENTITY, num_threads='all_cpus')\n        if self.original_img.count != 3:\n            self.layers = [rasterio.open(subd) for subd in self.original_img.subdatasets]\n                    \n        self.orig_size = self.original_img.shape\n\n        self.positions = self.get_positions()\n        \n        if rle is not None:\n            self.mask = enc2mask(rle, (self.orig_size[1], self.orig_size[0])) > 0\n        else:\n            self.mask = None\n\n    def __getitem__(self, idx):\n        \n        # Window\n        pos_x, pos_y = self.positions[idx]\n        x1, x2 = pos_x[0], pos_x[1]\n        y1, y2 = pos_y[0], pos_y[1]\n        window = Window.from_slices((x1, x2), (y1, y2))\n\n        # Retrieve slice\n        if self.original_img.count == 3:  # normal\n            img = self.original_img.read([1, 2, 3], window=window)\n            img = np.moveaxis(img, 0, -1)\n        else:  # with subdatasets/layers\n            img = np.zeros((self.tile_size, self.tile_size, 3), dtype=np.uint8)\n            for fl in range(3):\n                img[:, :, fl] = self.layers[fl].read(window=window) \n\n        # Downscale to tile size\n        img = cv2.resize(\n            img, (self.raw_tile_size, self.raw_tile_size), interpolation=cv2.INTER_AREA\n        )\n        img = self.transforms(image=img)[\"image\"]\n        \n        pos = np.array([pos_x[0], pos_x[1], pos_y[0], pos_y[1]])\n\n        return img, pos","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model","metadata":{}},{"cell_type":"markdown","source":"### Load weights","metadata":{}},{"cell_type":"code","source":"def load_model_weights(model, filename, verbose=1, cp_folder=\"\"):\n    \"\"\"\n    Loads the weights of a PyTorch model. The exception handles cpu/gpu incompatibilities.\n\n    Args:\n        model (torch model): Model to load the weights to.\n        filename (str): Name of the checkpoint.\n        verbose (int, optional): Whether to display infos. Defaults to 1.\n        cp_folder (str, optional): Folder to load from. Defaults to \"\".\n\n    Returns:\n        torch model: Model with loaded weights.\n    \"\"\"\n\n    if verbose:\n        print(f\"\\n -> Loading weights from {os.path.join(cp_folder,filename)}\\n\")\n    try:\n        model.load_state_dict(os.path.join(cp_folder, filename), strict=True)\n    except BaseException:\n        model.load_state_dict(\n            torch.load(os.path.join(cp_folder, filename), map_location=\"cpu\"),\n            strict=True,\n        )\n    return model","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Definition","metadata":{}},{"cell_type":"code","source":"import segmentation_models_pytorch\nfrom segmentation_models_pytorch.encoders import encoders\n\n\nDECODERS = [\"Unet\", \"Linknet\", \"FPN\", \"PSPNet\", \"DeepLabV3\", \"DeepLabV3Plus\", \"PAN\"]\nENCODERS = list(encoders.keys())\n\n\ndef define_model(\n    decoder_name, encoder_name, num_classes=1, activation=None, encoder_weights=\"imagenet\"\n):\n    \"\"\"\n    Loads a segmentation architecture\n\n    Args:\n        decoder_name (str): Decoder name.\n        encoder_name (str): Encoder name.\n        num_classes (int, optional): Number of classes. Defaults to 1.\n        pretrained : pretrained original weights\n    Returns:\n        torch model -- Pretrained model.\n    \"\"\"\n    assert decoder_name in DECODERS, \"Decoder name not supported\"\n    assert encoder_name in ENCODERS, \"Encoder name not supported\"\n\n    decoder = getattr(segmentation_models_pytorch, decoder_name)\n\n    model = decoder(\n        encoder_name,\n        encoder_weights=encoder_weights,\n        classes=num_classes,\n        activation=activation,\n    )\n    model.num_classes = num_classes\n\n    return model\n","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define & load","metadata":{}},{"cell_type":"code","source":"class Config:\n    def __init__(self, **entries):\n        self.__dict__.update(entries)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_models(cp_folder):\n    config = json.load(open(cp_folder + 'config.json', 'r'))\n    config = Config(**config)\n    \n    weights = sorted(glob.glob(cp_folder + \"*.pt\"))\n    models = []\n    \n    for weight in weights:\n        model = define_model(\n            config.decoder,\n            config.encoder,\n            num_classes=config.num_classes,\n            encoder_weights=None,\n        )\n        \n        model = load_model_weights(model, weight).to(DEVICE)\n        model.eval()\n        models.append(model)\n        \n#         break\n        \n    return models","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"markdown","source":"### Metric","metadata":{}},{"cell_type":"code","source":"def dice_scores_img(pred, truth, eps=1e-8):\n    \"\"\"\n    Dice metric for a single image as array.\n\n    Args:\n        pred (np array): Predictions.\n        truth (np array): Ground truths.\n        eps (float, optional): epsilon to avoid dividing by 0. Defaults to 1e-8.\n\n    Returns:\n        np array : dice value for each class\n    \"\"\"\n    pred = pred.reshape(-1) > 0\n    truth = truth.reshape(-1) > 0\n    intersect = (pred & truth).sum(-1)\n    union = pred.sum(-1) + truth.sum(-1)\n\n    dice = (2.0 * intersect + eps) / (union + eps)\n    return dice","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Tile weighting","metadata":{}},{"cell_type":"code","source":"def get_tile_weighting(size, sigma=1, alpha=1, eps=1e-6):\n    half = size // 2\n    w = np.ones((size, size), np.float32)\n\n    x = np.concatenate([np.mgrid[-half:0], np.mgrid[1: half + 1]])[:, None]\n    x = np.tile(x, (1, size))\n    x = half + 1 - np.abs(x)\n    y = x.T\n\n    w = np.minimum(x, y)\n    w = (w / w.max()) ** sigma\n    w = np.minimum(w, 1)\n\n    w = (w - np.min(w) + eps) / (np.max(w) - np.min(w) + eps)\n\n    w = np.where(w > alpha, 1, w)\n    w = w / alpha\n    w = np.clip(w, 1e-3, 1)\n\n    w = np.round(w, 3)\n    return w.astype(np.float16)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Predict","metadata":{}},{"cell_type":"code","source":"def predict_entire_mask(dataset, models, batch_size=32, tta=False):\n    loader = DataLoader(dataset, batch_size=batch_size, shuffle=False, pin_memory=True)\n\n    weighting = torch.from_numpy(get_tile_weighting(dataset.tile_size, sigma=1, alpha=1))\n    weighting_cuda = weighting.clone().cuda().unsqueeze(0)\n    weighting = weighting.cuda().half()\n\n    global_pred = torch.zeros(\n        (dataset.orig_size[0], dataset.orig_size[1]),\n        dtype=torch.half, device=\"cuda\"\n    )\n    global_counter = torch.zeros(\n        (dataset.orig_size[0], dataset.orig_size[1]),\n        dtype=torch.half, device=\"cuda\"\n    )\n\n    with torch.no_grad():\n        for img, pos in tqdm(loader):\n            img = img.to(\"cuda\")\n            _, _, h, w = img.shape\n            \n            model_preds = []\n            for model in models:\n                if model.num_classes == 1:\n                    pred = model(img).view(1, -1, h, w).sigmoid().detach()\n                else:\n                    pred = model(img)[:, 0].view(1, -1, h, w).sigmoid().detach()\n\n                if tta:\n                    for f in FLIPS:\n                        pred_flip = model(torch.flip(img, f))\n                        if model.num_classes == 2:\n                            pred_flip = pred_flip[:, :1]\n\n                        pred_flip = torch.flip(pred_flip, f).view(1, -1, h, w).sigmoid().detach()\n                        pred += pred_flip\n                    pred = torch.div(pred, len(FLIPS) + 1)\n\n                model_preds.append(pred)\n\n            pred = torch.cat(model_preds, 0).mean(0)\n\n            pred = torch.nn.functional.interpolate(\n                pred.unsqueeze(1), (dataset.tile_size, dataset.tile_size), mode='area'\n            ).squeeze(1)\n            \n            pred = (pred * weighting_cuda).half()\n\n            for tile_idx, (x0, x1, y0, y1) in enumerate(pos):\n                global_pred[x0: x1, y0: y1] += pred[tile_idx]\n                global_counter[x0: x1, y0: y1] += weighting\n\n    for i in range(len(global_pred)):\n        global_pred[i] = torch.div(global_pred[i], global_counter[i])\n\n    return global_pred","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Main","metadata":{}},{"cell_type":"code","source":"THRESHOLD = 0.5\n\nUSE_TTA = True # not DEBUG\nOVERLAP_FACTOR = 1.5\n\nCP_FOLDERS = [\n    \"../input/hubmap-cp/b1_2cfix/b1_2cfix/\",\n    \"../input/hubmap-cp/b1_last/b1_last/\"\n]\nCP_FOLDER = CP_FOLDERS[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(DATA_PATH + 'sample_submission.csv')\ndf_info = pd.read_csv(DATA_PATH + \"HuBMAP-20-dataset_information.csv\")\nrles = pd.read_csv(DATA_PATH + 'train.csv')\n\nconfig = json.load(open(CP_FOLDER + 'config.json', 'r'))\nconfig = Config(**config)\nconfig.overlap_factor = OVERLAP_FACTOR","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models = []\n\nfor cp_folder in CP_FOLDERS:\n    models += load_models(cp_folder)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for img in df['id'].unique():\n    if DEBUG:  # Check performances on a validation image\n        img = \"2f6ecfcdf\"  # check repro\n#         img = \"4ef6695ce\" # biggest img\n        IMG_PATH = DATA_PATH + \"train\"\n        models = [models[0], models[5]]\n                  \n    \n    print(f'\\n\\t Image {img}')\n    \n    if img == \"d488c759a\":\n        print('\\n - Using precomputed rle')\n        local_file_fc = '../input/hubmap-fast-submission/submission_0933_fc.csv'\n        df_local_fc = pd.read_csv(local_file_fc, index_col='id')\n        rle = df_local_fc.loc['d488c759a', 'predicted']\n        df.loc[df.id == img, 'predicted'] = rle\n\n        continue\n    \n    print(f'\\n - Building dataset')\n    \n    rle_truth = rles[rles['id'] == img][\"encoding\"] if DEBUG else None\n    \n    predict_dataset = InferenceEfficientDataset(\n        f\"{IMG_PATH}/{img}.tiff\",\n        rle=rle_truth,\n        overlap_factor=config.overlap_factor,\n        reduce_factor=config.reduce_factor,\n        tile_size=config.tile_size,\n        transforms=HE_preprocess(augment=False, visualize=False),\n    )\n    \n    print(f'\\n - Predicting masks')\n\n    global_pred = predict_entire_mask(\n        predict_dataset, models, batch_size=config.val_bs, tta=USE_TTA\n    )\n    \n    del predict_dataset\n    torch.cuda.empty_cache()\n    gc.collect()\n    \n    print('\\n - Encoding')\n    \n    global_pred_np = np.zeros(global_pred.size(), dtype=np.uint8)\n    for i in range(global_pred_np.shape[0]):\n        global_pred_np[i] = (global_pred[i] > THRESHOLD).cpu().numpy().astype(np.uint8)\n\n    rle = rle_encode_less_memory(global_pred_np)\n    df.loc[df.id == img, 'predicted'] = rle\n    \n    if DEBUG:\n        shape = df_info[df_info.image_file == img + \".tiff\"][['width_pixels', 'height_pixels']].values.astype(int)[0]\n        mask_truth = enc2mask(rle_truth, shape)\n        score = dice_scores_img(global_pred_np, mask_truth)\n        print(f\" -> Scored {score :.4f} with threshold {THRESHOLD:.2f}\")\n        break\n        \n    del global_pred, global_pred_np\n    torch.cuda.empty_cache()\n    gc.collect()\n\ndf.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}