{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":61446,"databundleVersionId":6962461,"sourceType":"competition"},{"sourceId":7158841,"sourceType":"datasetVersion","datasetId":4134545}],"dockerImageVersionId":30616,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Fast score computation using PyTorch\n\nI wrote a fast and memory-efficient code for the surface Dice score  *assuming tolerance 0*. The assumption greatly simplifies the code because there is no need to compute distances between surface elements, since distant surface elements do not contribute to the score. \n\nThe RAM usage is also significantly decreased by the order of 1/number of slices, because there are no need to compute pairs of distant surface elements, and therefore, only 2 slices of the 3D mask need to be on the RAM.\n\nThe computational complexity is O(W×H×D), where W, H, D are width, hight, and depth (number of slices) of the volume, respectively. RAM requirement is O(W×H) because only 2 slices are necessary simultaneously. These do not depend on the content of the segmentation at all; how many positive values or how fragmented they are (after rle are decoded to 2D array).\n\nI implemented with PyTorch. As a bonus, GPU can be used.\n\nThe scoreing of kidney_1_dense (1303×912×2279) took 20 minutes using the official code on my local machine, but reduced to 31 seconds with CPU and 3.1 seconds with GPU [3 minutes (CPU) and 16 seconds (GPU) with Kaggle notebook].\n\nThe score agree with the official implementation within float32 (~1e-7). Please let me know if you experience a discprancy between the scores or realise something wrong. Since I have just started this competition and just finished writing this code, there must be some bug and misunderstanding (hope they are not crutial).\n\n### Version 3:\nPyTorch version back to 1.90 can be used.","metadata":{}},{"cell_type":"markdown","source":"## 1. Score computation code\n\nThe code computes for one volume, like kidney_1_dense. If you want to compute for more than one volumes, like kidney_1 and kideny_3, compute separately and take the mean.\n\nSee the original paper about the surface Dice metric\n\n- https://arxiv.org/abs/1809.04430\n\nWhen you think about the surface Dice, imagine mask values 0/1 on 3-dimensional grid points, not in the cells (cubes or volume). Each cube in the grid has 8 corners and those corner points have the segmentation values 0 or 1. There are 256 patterns of 0/1 for 8 corners, and for each pattern they define a 2-dimmensional surface in the cube, known as the marching cubes.\n\nhttps://en.wikipedia.org/wiki/Marching_cubes\n\nTrue-positiveness is same as the ordinary Dice score, If both prediction and ground truth have a non-zero surface in the same cube, that cube is true positive. The 8-corner pattern does not have to be the same to be true positive.\n\nThe difference between the volume-Dice and surface-Dice is that the surface area is zero if that cube is completely inside the volume (i.e., all 8 corners have value 1). Only the boundaries contribute to the surface Dice. Also, not all boundary cubes have the same weights, the weight is the marching-cube surface area.\n\n$$\n\\mathtt{surface\\_dice}= \\frac{\\sum_{tp} S_\\mathrm{pred} + \\sum_{tp} S_\\mathrm{true}}{\\sum_\\mathrm{all} S_\\mathrm{pred} + \\sum_\\mathrm{all} S_\\mathrm{true}}\n$$\n\nThe sum is over all cubes in the volume and S is the surface area in those cubes. $\\sum_\\mathrm{tp}$ is a sum over all true positive cubes, where both prediction and ground truth have non-zero areas.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport sys\n\nimport time\nimport torch\nimport torch.nn as nn\nfrom PIL import Image\n\nsys.path.append('/kaggle/input/sennet-score/src')\nfrom official_metric import create_table_neighbour_code_to_surface_area\n\ndi = '/kaggle/input/blood-vessel-segmentation'\ndevice = torch.device('cuda')  # can be 'cpu'\n\n# PyTorch version dependence on index data type\ntorch_ver_major = int(torch.__version__.split('.')[0])\ndtype_index = torch.int32 if torch_ver_major >= 2 else torch.long","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:07:36.286531Z","iopub.execute_input":"2023-12-13T00:07:36.287439Z","iopub.status.idle":"2023-12-13T00:07:40.772285Z","shell.execute_reply.started":"2023-12-13T00:07:36.287403Z","shell.execute_reply":"2023-12-13T00:07:40.771309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def rle_decode(mask_rle: str, shape: tuple) -> np.array:\n    \"\"\"\n    Decode rle string\n    https://www.kaggle.com/code/paulorzp/run-length-encode-and-decode/script\n    https://www.kaggle.com/stainsby/fast-tested-rle\n\n    Args:\n      mask_rle: run length (rle) as string\n      shape: (height, width) of the mask\n\n    Returns:\n      array[uint8], 1 - mask, 0 - background\n    \"\"\"\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0] * shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n    return img.reshape(shape)\n\n\ndef compute_area(y: list, unfold: nn.Unfold, area: torch.Tensor) -> torch.Tensor:\n    \"\"\"\n    Args:\n      y (list[Tensor]): A pair of consecutive slices of mask\n      unfold: nn.Unfold(kernel_size=(2, 2), padding=1)\n      area (Tensor): surface area for 256 patterns (256, )\n\n    Returns:\n      Surface area of surface in 2x2x2 cube\n    \"\"\"\n    # Two layers of segmentation masks\n    yy = torch.stack(y, dim=0).to(torch.float16).unsqueeze(0)\n    # (batch_size=1, nch=2, H, W) \n    # bit (0/1) but unfold requires float\n\n    # unfold slides through the volume like a convolution\n    # 2x2 kernel returns 8 values (2 channels * 2x2)\n    cubes_float = unfold(yy).squeeze(0)  # (8, n_cubes)\n\n    # Each of the 8 values are either 0 or 1\n    # Convert those 8 bits to one uint8\n    cubes_byte = torch.zeros(cubes_float.size(1), dtype=dtype_index, device=device)\n    # indices are required to be int32 or long for area[cube_byte] below, not uint8\n    # Can be int32 for torch 2.0.0, int32 raise IndexError in torch 1.13.1.\n    \n    for k in range(8):\n        cubes_byte += cubes_float[k, :].to(dtype_index) << k\n\n    # Use area lookup table: pattern index -> area [float]\n    cubes_area = area[cubes_byte]\n\n    return cubes_area\n\n\ndef compute_surface_dice_score(submit: pd.DataFrame, label: pd.DataFrame) -> float:\n    \"\"\"\n    Compute surface Dice score for one 3D volume\n\n    submit (pd.DataFrame): submission file with id and rle\n    label (pd.DataFrame): ground truth id, rle, and also image height, width\n    \"\"\"\n    # submit and label must contain exact same id in same order\n    assert (submit['id'] == label['id']).all()\n    assert len(label) > 0\n\n    # All height, width must be the same\n    len(label['height'].unique()) == 1\n    len(label['width'].unique()) == 1\n\n    # Surface area lookup table: Tensor[float32] (256, )\n    area = create_table_neighbour_code_to_surface_area((1, 1, 1))\n    area = torch.from_numpy(area).to(device)  # torch.float32\n\n    # Slide through the volume like a convolution\n    unfold = torch.nn.Unfold(kernel_size=(2, 2), padding=1)\n\n    r = label.iloc[0]\n    h, w = r['height'], r['width']\n    n_slices = len(label)\n\n    # Padding before first slice\n    y0 = y0_pred = torch.zeros((h, w), dtype=torch.uint8, device=device)\n\n    num = 0     # numerator of surface Dice\n    denom = 0   # denominator\n    for i in range(n_slices + 1):\n        # Load one slice\n        if i < n_slices:\n            r = label.iloc[i]\n            y1 = rle_decode(r['rle'], (h, w))\n            y1 = torch.from_numpy(y1).to(device)\n\n            r = submit.iloc[i]\n            y1_pred = rle_decode(r['rle'], (h, w))\n            y1_pred = torch.from_numpy(y1_pred).to(device)\n        else:\n            # Padding after the last slice\n            y1 = y1_pred = torch.zeros((h, w), dtype=torch.uint8, device=device)\n\n        # Compute the surface area between two slices (n_cubes,)\n        area_pred = compute_area([y0_pred, y1_pred], unfold, area)\n        area_true = compute_area([y0, y1], unfold, area)\n\n        # True positive cube indices\n        idx = torch.logical_and(area_pred > 0, area_true > 0)\n\n        # Surface dice numerator and denominator\n        num += area_pred[idx].sum() + area_true[idx].sum()\n        denom += area_pred.sum() + area_true.sum()\n\n        # Next slice\n        y0 = y1\n        y0_pred = y1_pred\n\n    dice = num / denom.clamp(min=1e-8)\n    return dice.item()","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:07:40.774392Z","iopub.execute_input":"2023-12-13T00:07:40.774905Z","iopub.status.idle":"2023-12-13T00:07:40.794260Z","shell.execute_reply.started":"2023-12-13T00:07:40.774871Z","shell.execute_reply":"2023-12-13T00:07:40.793238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. Use the scoring function\n\n### Prediction\n\n`submission1.csv` is my prediction for kidney_1_dense using trivial 2D U-Net similar to what [yumeneko has shared](https://www.kaggle.com/code/kashiwaba/sennet-hoa-train-unet-simple-baseline). The score is irrelevant because this is trained with kidney_1 and applied to kidney_1 (leakage).","metadata":{}},{"cell_type":"code","source":"submit = pd.read_csv('/kaggle/input/sennet-score/submission1.csv')\nsubmit.head(n=2)","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:07:40.795249Z","iopub.execute_input":"2023-12-13T00:07:40.795602Z","iopub.status.idle":"2023-12-13T00:07:41.246380Z","shell.execute_reply.started":"2023-12-13T00:07:40.795570Z","shell.execute_reply":"2023-12-13T00:07:41.245444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Need height and weight columns; this is same for the official code.","metadata":{}},{"cell_type":"code","source":"def add_size_columns(df: pd.DataFrame):\n    \"\"\"\n    df (DataFrame): including id column, e.g., kidney_1_dense_0000\n    \"\"\"\n    widths = []\n    heights = []\n    subdirs = []\n    nums = []\n    for i, r in df.iterrows():\n        file_id = r['id']\n        subdir = file_id[:-5]    # kidney_1_dense\n        file_num = file_id[-4:]  # 0000\n\n        filename = '%s/train/%s/images/%s.tif' % (di, subdir, file_num)\n        img = Image.open(filename)\n        w, h = img.size\n        widths.append(w)\n        heights.append(h)\n        subdirs.append(subdir)\n        nums.append(file_num)\n\n    df['width'] = widths\n    df['height'] = heights\n    df['image_id'] = subdirs\n    df['slice_id'] = nums","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:07:41.249165Z","iopub.execute_input":"2023-12-13T00:07:41.249530Z","iopub.status.idle":"2023-12-13T00:07:41.256769Z","shell.execute_reply.started":"2023-12-13T00:07:41.249497Z","shell.execute_reply":"2023-12-13T00:07:41.255895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Label for kidney_1_dense only\nimage_id = 'kidney_1_dense'\nlabel = pd.read_csv(di + '/train_rles.csv')\nidx = label['id'].str.startswith(image_id)\nlabel = label[idx]\nassert len(label) > 0\n\n# Add height, width columns\nadd_size_columns(label)\n\nlabel.head(n=2)","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:07:41.257871Z","iopub.execute_input":"2023-12-13T00:07:41.258202Z","iopub.status.idle":"2023-12-13T00:08:16.369562Z","shell.execute_reply.started":"2023-12-13T00:07:41.258176Z","shell.execute_reply":"2023-12-13T00:08:16.368614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Compute surface Dice score\nscore = compute_surface_dice_score(submit, label)\nscore","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:08:16.370709Z","iopub.execute_input":"2023-12-13T00:08:16.370985Z","iopub.status.idle":"2023-12-13T00:08:33.311398Z","shell.execute_reply.started":"2023-12-13T00:08:16.370961Z","shell.execute_reply":"2023-12-13T00:08:33.310500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Compare with the official implementation\n\nofficial_metric.py is a copy of the notebook (Version 24: 6 Dec 2023):\n\nhttps://www.kaggle.com/code/metric/surface-dice-metric/notebook\n\nDue to insufficient memory and my patience, I reduce the number of slices to 100.","metadata":{}},{"cell_type":"code","source":"begin, end = 1000, 1100\nlabel = label.iloc[begin:end]\nsubmit = submit.iloc[begin:end]","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:08:33.312782Z","iopub.execute_input":"2023-12-13T00:08:33.313169Z","iopub.status.idle":"2023-12-13T00:08:33.318700Z","shell.execute_reply.started":"2023-12-13T00:08:33.313134Z","shell.execute_reply":"2023-12-13T00:08:33.317795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Recompute for the subvolume\nscore = compute_surface_dice_score(submit, label)\nscore","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:08:33.320066Z","iopub.execute_input":"2023-12-13T00:08:33.320374Z","iopub.status.idle":"2023-12-13T00:08:34.073367Z","shell.execute_reply.started":"2023-12-13T00:08:33.320349Z","shell.execute_reply":"2023-12-13T00:08:34.072366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Again, this is a score with leakage (or train score).","metadata":{}},{"cell_type":"code","source":"import official_metric","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:08:34.074642Z","iopub.execute_input":"2023-12-13T00:08:34.074951Z","iopub.status.idle":"2023-12-13T00:08:34.078966Z","shell.execute_reply.started":"2023-12-13T00:08:34.074926Z","shell.execute_reply":"2023-12-13T00:08:34.078157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This will take 2 minutes for 100 slices.","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:32:55.602157Z","iopub.execute_input":"2023-12-09T06:32:55.602881Z","iopub.status.idle":"2023-12-09T06:32:55.614258Z","shell.execute_reply.started":"2023-12-09T06:32:55.602844Z","shell.execute_reply":"2023-12-09T06:32:55.612303Z"}}},{"cell_type":"code","source":"%%time\nscore_official = official_metric.score(label, submit,\n                                       row_id_column_name='id',\n                                       rle_column_name='rle',\n                                       tolerance=0.0,\n                                       image_id_column_name='image_id',\n                                       slice_id_column_name='slice_id')\nscore_official","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:08:34.081904Z","iopub.execute_input":"2023-12-13T00:08:34.082598Z","iopub.status.idle":"2023-12-13T00:10:09.591737Z","shell.execute_reply.started":"2023-12-13T00:08:34.082561Z","shell.execute_reply":"2023-12-13T00:10:09.590822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score - score_official","metadata":{"execution":{"iopub.status.busy":"2023-12-13T00:10:09.592860Z","iopub.execute_input":"2023-12-13T00:10:09.593181Z","iopub.status.idle":"2023-12-13T00:10:09.599002Z","shell.execute_reply.started":"2023-12-13T00:10:09.593157Z","shell.execute_reply":"2023-12-13T00:10:09.598122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Acceptable ~1e-7 difference for float32.","metadata":{}}]}