{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":61446,"databundleVersionId":6962461,"sourceType":"competition"},{"sourceId":162928788,"sourceType":"kernelVersion"}],"dockerImageVersionId":30627,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<style>\n    .header {font-size: 40px}\n</style>\n<p style=\"font-size: 40px\">\n    <img align=\"right\" width=\"400\" height=\"200\" src=\"https://www.kaggle.com/competitions/61446/images/header\">\n    <b>SenNet + HOA - Hacking the Human Vasculature in 3D</b>\n    <p style=\"font-size: 18px\">Segment vasculature in 3D scans of human kidney</p>\n</p>","metadata":{}},{"cell_type":"markdown","source":"This solution was implemented as part of a blood vessel segmentation competition organized by the Common Fund’s Cellular Senescence Network (SenNet) Programm in cooperation with the Human Organ Atlas (HOA). \n - Competition overview page: [SenNet + HOA - Hacking the Human Vasculature in 3D](https://www.kaggle.com/competitions/blood-vessel-segmentation)\n - Competition dataset is [here](https://www.kaggle.com/competitions/blood-vessel-segmentation/data)","metadata":{}},{"cell_type":"markdown","source":"**NB**: This notebook contains only the inference stage, the training stage is [here](https://www.kaggle.com/code/pib73nl/sennet-hoa-bvs-12th-place-solution-train)","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport os, glob\nimport matplotlib.pyplot as plt\nimport gc\nfrom tqdm import tqdm\nimport tensorflow as tf\nimport skimage as ski\nfrom skimage import io\n\nprint(f'TF version - {tf.__version__}')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-15T06:51:45.762077Z","iopub.execute_input":"2024-02-15T06:51:45.762733Z","iopub.status.idle":"2024-02-15T06:51:58.951341Z","shell.execute_reply.started":"2024-02-15T06:51:45.762695Z","shell.execute_reply":"2024-02-15T06:51:58.950379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# detect hardware, return appropriate distribution strategy\n# for this competition use GPU, TPU is under restrictions\ntry:\n    tpu = tf.distribute.cluster_resolver.TPUClusterResolver()\n    tf.config.experimental_connect_to_cluster(tpu)\n    tf.tpu.experimental.initialize_tpu_system(tpu)\n    strategy = tf.distribute.TPUStrategy(tpu)\n    print('Running on TPU ', tpu.master())\nexcept ValueError:\n    strategy = tf.distribute.MirroredStrategy() # for CPU/GPU or multi-GPU machines\n\nprint(\"Number of replicas: \", strategy.num_replicas_in_sync)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:51:58.953095Z","iopub.execute_input":"2024-02-15T06:51:58.953592Z","iopub.status.idle":"2024-02-15T06:52:00.356758Z","shell.execute_reply.started":"2024-02-15T06:51:58.953564Z","shell.execute_reply":"2024-02-15T06:52:00.355776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MAIN_PATH = '/kaggle/input/blood-vessel-segmentation'\nTEST_PATH = MAIN_PATH + '/test'\nCHECK_PATH = MAIN_PATH + '/train/kidney_1_dense'\nWEIGHTS_PATH = '/kaggle/input/sennet-hoa-bvs-12th-place-solution-train/tmp/checkpoint'\n\nINPUT_SHAPE = [1024, 704]\n# INPUT_SHAPE = [896, 576]\n# INPUT_SHAPE = [768, 512]\n# INPUT_SHAPE = [576, 384]\n# INPUT_SHAPE = [384, 256]\n\nDEBUG = True","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.358035Z","iopub.execute_input":"2024-02-15T06:52:00.358341Z","iopub.status.idle":"2024-02-15T06:52:00.363432Z","shell.execute_reply.started":"2024-02-15T06:52:00.358314Z","shell.execute_reply":"2024-02-15T06:52:00.362351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Building a model and loading weights","metadata":{}},{"cell_type":"markdown","source":"Re-building a model and loading weights is used, since loading the saved model failed due to errors (different versions of the TF or something else...)","metadata":{}},{"cell_type":"code","source":"# Encoder Utilities\n\ndef conv2d_block(input_tensor, n_filters, kernel_size = 3):\n    \n    '''\n    Adds 2 convolutional layers with the parameters passed to it\n\n    Args:\n    input_tensor (tensor) -- the input tensor\n    n_filters (int) -- number of filters\n    kernel_size (int) -- kernel size for the convolution\n\n    Returns:\n    tensor of output features\n    '''\n    # first layer\n    x = input_tensor\n    for _ in range(2):\n        x = tf.keras.layers.Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size), \n                                   kernel_initializer = 'random_uniform', padding = 'same')(x)\n        x = tf.keras.layers.Activation('relu')(x)\n\n    return x\n\n\ndef encoder_block(inputs, n_filters=64, pool_type = 'max', pool_size=(2,2), dropout=0.5, normalize = True):\n    '''\n    Adds two convolutional blocks and then perform down sampling on output of convolutions.\n\n    Args:\n    input_tensor (tensor) -- the input tensor\n    n_filters (int) -- number of filters\n    pool_type (str) -- type of pooling ('max', 'avg')\n    pool_size (tuple) -- size of a pooling window\n    dropout (float) -- dropout rate\n    normalize (bool) -- whether to use normalization\n\n    Returns:\n    f - the output features of the convolution block \n    p - the pooled features with dropout and normalization\n    '''\n\n    f = conv2d_block(inputs, n_filters=n_filters)\n    \n    if pool_type == 'max':\n        p = tf.keras.layers.MaxPooling2D(pool_size=pool_size)(f)\n    elif pool_type == 'avg':\n        p = tf.keras.layers.AveragePooling2D(pool_size=pool_size)(f)\n    \n    if normalize:\n        p = tf.keras.layers.BatchNormalization()(p)\n        \n    p = tf.keras.layers.Dropout(dropout)(p)\n\n    return f, p\n\n\ndef encoder(inputs):\n    '''\n    This function defines the encoder or downsampling path.\n\n    Args:\n    inputs (tensor) -- batch of input images\n\n    Returns:\n    p4 - the output maxpooled features of the last encoder block\n    (f1, f2, f3, f4) - the output features of all the encoder blocks\n    '''\n    f1, p1 = encoder_block(inputs, n_filters=64, pool_type='avg')\n    f2, p2 = encoder_block(p1, n_filters=128, pool_type='avg')\n    f3, p3 = encoder_block(p2, n_filters=256, pool_type='avg', normalize=False)\n    f4, p4 = encoder_block(p3, n_filters=512, pool_type='avg', normalize=False)\n\n    return p4, (f1, f2, f3, f4)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.366193Z","iopub.execute_input":"2024-02-15T06:52:00.366464Z","iopub.status.idle":"2024-02-15T06:52:00.385132Z","shell.execute_reply.started":"2024-02-15T06:52:00.366439Z","shell.execute_reply":"2024-02-15T06:52:00.383878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def bottleneck(inputs):\n    '''\n    This function defines the bottleneck convolutions to extract more features before the upsampling layers.\n    '''\n\n    bottle_neck = conv2d_block(inputs, n_filters=1024)\n\n    return bottle_neck","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.386473Z","iopub.execute_input":"2024-02-15T06:52:00.386814Z","iopub.status.idle":"2024-02-15T06:52:00.399777Z","shell.execute_reply.started":"2024-02-15T06:52:00.386773Z","shell.execute_reply":"2024-02-15T06:52:00.398893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Decoder Utilities\n\ndef decoder_block(inputs, conv_output, n_filters=64, kernel_size=3, strides=(2,2), dropout=0.5, normalize = True):\n    '''\n    defines the one decoder block of the UNet\n\n    Args:\n    inputs (tensor) -- batch of input features\n    conv_output (tensor) -- features from an encoder block\n    n_filters (int) -- number of filters\n    kernel_size (int) -- kernel size\n    strides (int) -- strides for the deconvolution/upsampling\n    dropout (float) -- dropout rate\n    normalize (bool) -- whether to use normalization\n\n    Returns:\n    c (tensor) -- output features of the decoder block\n    '''\n    u = tf.keras.layers.Conv2DTranspose(n_filters, kernel_size, strides = strides, padding = 'same')(inputs)\n    c = tf.keras.layers.concatenate([u, conv_output])\n    \n    if normalize:\n        c = tf.keras.layers.BatchNormalization()(c)\n    \n    c = tf.keras.layers.Dropout(dropout)(c)\n    \n    c = conv2d_block(c, n_filters, kernel_size=kernel_size)\n\n    return c\n\n\ndef decoder(inputs, convs, output_channels):\n    '''\n    Defines the decoder of the UNet chaining together several decoder blocks. \n\n    Args:\n    inputs (tensor) -- batch of input features\n    convs (tuple) -- features from the encoder blocks\n    output_channels (int) -- number of classes in the label map\n\n    Returns:\n    outputs (tensor) -- the pixel wise label map of the image\n    '''\n\n    f1, f2, f3, f4 = convs\n\n    c6 = decoder_block(inputs, f4, n_filters=512, normalize=False)\n    c7 = decoder_block(c6, f3, n_filters=256, normalize=False)\n    c8 = decoder_block(c7, f2, n_filters=128)\n    c9 = decoder_block(c8, f1, n_filters=64)\n\n    outputs = tf.keras.layers.Conv2D(output_channels, (1, 1), activation='sigmoid')(c9)\n\n    return outputs","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.400955Z","iopub.execute_input":"2024-02-15T06:52:00.401302Z","iopub.status.idle":"2024-02-15T06:52:00.414133Z","shell.execute_reply.started":"2024-02-15T06:52:00.401268Z","shell.execute_reply":"2024-02-15T06:52:00.413132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"OUTPUT_CHANNELS = 1\n\ndef unet():\n    '''\n    Defines the UNet by connecting the encoder, bottleneck and decoder.\n    '''\n\n    # specify the input shape\n    inputs = tf.keras.layers.Input(shape=(*INPUT_SHAPE, 1))\n\n    # feed the inputs to the encoder\n    encoder_output, convs = encoder(inputs)\n\n    # feed the encoder output to the bottleneck\n    bottle_neck = bottleneck(encoder_output)\n\n    # feed the bottleneck and encoder block outputs to the decoder\n    # specify the number of classes via the `output_channels` argument\n    outputs = decoder(bottle_neck, convs, output_channels=OUTPUT_CHANNELS)\n\n    # create the model\n    model = tf.keras.Model(inputs=inputs, outputs=outputs)\n\n    return model","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.415471Z","iopub.execute_input":"2024-02-15T06:52:00.416099Z","iopub.status.idle":"2024-02-15T06:52:00.428138Z","shell.execute_reply.started":"2024-02-15T06:52:00.416049Z","shell.execute_reply":"2024-02-15T06:52:00.426807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with strategy.scope():\n    model = unet()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:00.430057Z","iopub.execute_input":"2024-02-15T06:52:00.430877Z","iopub.status.idle":"2024-02-15T06:52:01.522029Z","shell.execute_reply.started":"2024-02-15T06:52:00.430842Z","shell.execute_reply":"2024-02-15T06:52:01.520986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# loading weights\nmodel.load_weights(WEIGHTS_PATH)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:01.523521Z","iopub.execute_input":"2024-02-15T06:52:01.523925Z","iopub.status.idle":"2024-02-15T06:52:03.460077Z","shell.execute_reply.started":"2024-02-15T06:52:01.523871Z","shell.execute_reply":"2024-02-15T06:52:03.458949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"code","source":"def read_image_label(im_path, lb_path = None):\n    \"\"\"\n    Pre-process images and labels (if exist) and remember metadata to restore the original size\n    \"\"\"\n    image = io.imread(im_path)\n    label = io.imread(lb_path) if lb_path else None\n    # for restoring original size and taking ID\n    meta = {\n            'orig_shape': image.shape,\n            'source':im_path.split('/')[-3::2], # last three elements [kidney_*, images, *.tif] with step 2 - [kidney_*, *.tif]\n           }\n    \n    # normalizing before applying roi\n    image = (image - image.mean()) / image.std()\n    image, label, pad_width = apply_roi(image, label)\n    # normalizing after applying roi\n    image = (image - image.mean()) / image.std()\n    \n    meta['roi_shape'] = image.shape\n    meta['pad_width'] = pad_width\n    image = ski.transform.resize(image, INPUT_SHAPE).astype(np.float16)\n\n    if isinstance(label, np.ndarray):\n        label = (label / 255).astype(np.uint8)\n        label = ski.transform.resize(label, INPUT_SHAPE, order=0) \n        # the following construction is not used further (but should be used); left here for further improvement\n        class_weights = [(1 - label.mean()).astype(np.float32), label.mean().astype(np.float32)] if label.mean() != 0 else [1.0, 0.0]\n    else:\n        class_weights = None\n    \n    \n    return image, label, class_weights, meta\n\ndef apply_roi(image, label=None):\n    \"\"\"\n    Exclusion of uninformative image fields \n    \"\"\"\n    # just throw out rows and columnt with low std\n    row_mask = image.std(axis=1)>0.22\n    clmn_mask = image.std(axis=0)>0.22\n    \n    # cleaning up the of noize of this approach and taking a solid region\n    row_mask = cleaning_mask(row_mask)\n    clmn_mask = cleaning_mask(clmn_mask)\n    \n    image = image[row_mask,:][:, clmn_mask]\n    label = label[row_mask,:][:, clmn_mask] if isinstance(label, np.ndarray) else None \n    \n    # remember the size of the pads for subsequent correct restoration\n    row_pad = (row_mask.argmax(), row_mask[::-1].argmax())\n    clmn_pad = (clmn_mask.argmax(), clmn_mask[::-1].argmax())\n    \n    return image, label, (row_pad, clmn_pad)\n\ndef cleaning_mask(mask):\n    \"\"\"\n    Selecting a solid region from a noisy mask\n    \"\"\"\n    # if frame starts from the first element or finishes at the last\n    mask[0] = False\n    mask[-1] = False\n    \n    # taking edges of frames\n    frames = np.nonzero(mask[:-1]!=mask[1:])[0]\n    # taking length of frames\n    delta = frames[1:]-frames[:-1]\n    # taking index of max len frame\n    max_solid_block_begin = np.argmax(delta)\n    # other is garbage\n    garbage = np.delete(frames, [max_solid_block_begin, max_solid_block_begin+1])\n    # clearing the mask\n    for a, b in zip(garbage[::2], garbage[1::2]):\n        mask[a+1:b+1] = False\n    \n    return mask","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:03.463743Z","iopub.execute_input":"2024-02-15T06:52:03.464154Z","iopub.status.idle":"2024-02-15T06:52:03.481959Z","shell.execute_reply.started":"2024-02-15T06:52:03.464115Z","shell.execute_reply":"2024-02-15T06:52:03.480921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resize2original_shape(image, roi_shape, pad_width):\n    \"\"\"\n    Restoring the label to its original size\n    \"\"\"\n    # returning size to shape after applying roi\n    resized2roi = ski.transform.resize(np.squeeze(image), roi_shape, order=0)\n    # adding zero pads to restore original size\n    resized2orig = np.pad(resized2roi, pad_width, mode='constant', constant_values=0)\n    \n    return resized2orig\n\ndef rle_encode(mask):\n    \"\"\"\n    rle encoder (thanks to the community)\n    \"\"\"\n    pixel = mask.flatten()\n    pixel[pixel<0.5] = 0\n    pixel[pixel>0.5] = 1 \n    pixel = np.concatenate([[0], pixel, [0]])\n    run = np.where(pixel[1:] != pixel[:-1])[0] + 1\n    run[1::2] -= run[::2]\n    rle = ' '.join(str(r) for r in run)\n    if rle == '':\n        rle = '1 1'\n    return rle\n\ndef prediction(image_path, labeled=False, resize=False, rle=False):\n    \"\"\"\n    Full cycle: read & predprocess images & labels (if exist);\n    prediction, restoring labels to its original size, rle encoding \n    \"\"\"\n    \n    predict = {}\n    \n    for im_path in tqdm(image_path, ncols=80, total=len(image_path)):\n        # reading & preprocessing\n        image, label, sample_weight, meta = read_image_label(im_path, im_path.replace('image','label') if labeled else None)\n\n        id = f\"{meta['source'][0]}_{os.path.splitext(meta['source'][1])[0]}\"\n        # prediction\n        pred = model.predict(image[np.newaxis, ...], verbose=0)\n        \n        # restoring to original size\n        if resize:\n            pred = resize2original_shape(np.squeeze(pred), meta['roi_shape'], meta['pad_width'])\n            if labeled:\n                label = (io.imread(im_path.replace('image','label')) / 255).astype(np.uint8) \n        \n        # rle encoding\n        if rle:\n            pred = rle_encode(pred)\n        \n        # if labeled images are given (for checking the solution)\n        if labeled:\n            predict[id] = {id: value for id, value in zip(['label', 'pred'], [label, pred])}\n        else:\n            predict[id] = pred\n        \n    return predict","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:03.483336Z","iopub.execute_input":"2024-02-15T06:52:03.483719Z","iopub.status.idle":"2024-02-15T06:52:03.501422Z","shell.execute_reply.started":"2024-02-15T06:52:03.483683Z","shell.execute_reply":"2024-02-15T06:52:03.500390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_path_tst = sorted(glob.glob(f'{TEST_PATH}/*/*/*.tif'))\n# for submition need to set resize=True, rle=True\npred_dict = prediction(image_path_tst, resize=True, rle=True)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:03.502728Z","iopub.execute_input":"2024-02-15T06:52:03.503093Z","iopub.status.idle":"2024-02-15T06:52:16.395458Z","shell.execute_reply.started":"2024-02-15T06:52:03.503058Z","shell.execute_reply":"2024-02-15T06:52:16.394441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame.from_dict(pred_dict, orient='index', columns=['rle'])\nsubmission.head(10)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:16.396848Z","iopub.execute_input":"2024-02-15T06:52:16.397264Z","iopub.status.idle":"2024-02-15T06:52:16.419808Z","shell.execute_reply.started":"2024-02-15T06:52:16.397224Z","shell.execute_reply":"2024-02-15T06:52:16.418960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('submission.csv', index_label='id', columns=['rle'])","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:16.420924Z","iopub.execute_input":"2024-02-15T06:52:16.421211Z","iopub.status.idle":"2024-02-15T06:52:16.434889Z","shell.execute_reply.started":"2024-02-15T06:52:16.421187Z","shell.execute_reply":"2024-02-15T06:52:16.434135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Checking solution","metadata":{}},{"cell_type":"code","source":"# let's take some labeled images for checking the solution\nif DEBUG:\n    image_path_chck = sorted(glob.glob(f'{CHECK_PATH}/images/*.tif'))\n    step = len(image_path_chck) // 5\n    image_path_chck = (\n        image_path_chck[0:3] + \n        image_path_chck[step:step+3] + \n        image_path_chck[step*2:step*2+3] + \n        image_path_chck[step*3:step*3+3] +\n        image_path_chck[step*4:step*4+3] +\n        image_path_chck[step*5:step*5+3]\n    ) \n\n    pred_dict = prediction(image_path_chck, labeled=True, resize=True)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:16.435888Z","iopub.execute_input":"2024-02-15T06:52:16.436194Z","iopub.status.idle":"2024-02-15T06:52:26.251317Z","shell.execute_reply.started":"2024-02-15T06:52:16.436168Z","shell.execute_reply":"2024-02-15T06:52:26.250334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if DEBUG:\n    plt.figure(figsize=(30, 300))\n\n    rows = len(pred_dict.keys())\n\n    for i, id in zip(range(1, rows*3, 3), pred_dict.keys()):\n\n        lbl = pred_dict[id]['label']\n        pred = np.squeeze(pred_dict[id]['pred'])\n        pred[pred<0.5] = 0\n        pred[pred>0.5] = 1\n        fp = ((lbl-pred)<0).sum()\n        fn = ((lbl-pred)>0).sum()\n        tp = ((lbl+pred)==2).sum()\n        tn = ((lbl+pred)==0).sum()\n\n        ax = plt.subplot(rows, 3, i)\n        plt.imshow(lbl, cmap='gray')\n        plt.axis(\"off\")\n        plt.title('ground truth', fontdict={'fontsize':25})\n\n        ax = plt.subplot(rows, 3, i + 1)\n        plt.imshow(pred, cmap='gray')\n        plt.axis(\"off\")\n        plt.title(f'prediction\\nfp - {fp}, fn - {fn}\\ntp - {tp}, tn - {tn}', fontdict={'fontsize':25})\n        \n        ax = plt.subplot(rows, 3, i+2)\n        plt.imshow(lbl, cmap='gray')\n        plt.imshow(pred, cmap='coolwarm', alpha=0.5)\n        plt.axis(\"off\")\n        plt.title('overlap', fontdict={'fontsize':25})","metadata":{"execution":{"iopub.status.busy":"2024-02-15T06:52:26.252642Z","iopub.execute_input":"2024-02-15T06:52:26.252956Z","iopub.status.idle":"2024-02-15T06:52:48.800463Z","shell.execute_reply.started":"2024-02-15T06:52:26.252918Z","shell.execute_reply":"2024-02-15T06:52:48.799139Z"},"trusted":true},"execution_count":null,"outputs":[]}]}