{"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":7183197,"sourceType":"datasetVersion","datasetId":4152203},{"sourceId":155342803,"sourceType":"kernelVersion"}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Post Competition notebook 2D to 3D UNet for Private LB 0.611960\n## with connected components 3d (cc3d) post processing\n\nDuring the competition there were discussions on post processing and connected components to obtain the best tree like structure for the human kidney vasculature. The package connected components 3d (cc3d) was  one idea which looked promising.  \n\nPost Competition testing to see how well this could perform was undertaken and this notebook has the best obtained results on a resnet50 2D to 3D UNet. It was found with cc3d.dust that a 6-connected neighborhood worked better than higher values 18 or 26. And best threshold 6552 to remove objects with fewer than the threshold voxels. (note: threshold 6474 - 6630 performed equally well.)  Prediction probability threshold of 0.29 was better than 0.3 or 0.4.   \n\nPerhaps even better results could be found with other models.  Private LB 0.611960 would have been high in the silver range in this competition. Using connected components 3d could be a useful strategy in future work. \n\n\nThe original source of code and model are from hengck23, please refer to: \n\nhttps://www.kaggle.com/code/hengck23/submission-example-for-2d-to-3d-unet\n\nhttps://www.kaggle.com/datasets/hengck23/blood-vessel-segmentation-01\n\nand competition discussion thread on experiments and insights:\n\nhttps://www.kaggle.com/competitions/blood-vessel-segmentation/discussion/456118\n\nPackage link:\n\nhttps://pypi.org/project/connected-components-3d/\n\n\"cc3d is an implementation of connected components in three dimensions using a 26, 18, or 6-connected neighborhood in 3D or 4 and 8-connected in 2D. This package uses a 3D variant of the two pass method by Rosenfeld and Pflatz augmented with Union-Find and a decision tree based on the 2D 8-connected work of Wu, Otoo, and Suzuki. This implementation is compatible with images containing many different labels, not just binary images. It also supports continuously valued images such as grayscale microscope images with an algorithm that joins together nearby values.\"\n\ncc3d.dust refers to  \n\"Removal of small objects (\"dust\")\"","metadata":{}},{"cell_type":"markdown","source":"# Package Install","metadata":{}},{"cell_type":"code","source":"!mkdir -p /tmp/pip/cache/\n!cp /kaggle/input/sennet-hoa-pkgs/connected_components_3d-3.12.4-cp310-cp310-linux_x86_64.whl /tmp/pip/cache/\n!ls /tmp/pip/cache/","metadata":{"execution":{"iopub.status.busy":"2023-12-17T08:09:35.951399Z","iopub.execute_input":"2023-12-17T08:09:35.951813Z","iopub.status.idle":"2023-12-17T08:09:39.405591Z","shell.execute_reply.started":"2023-12-17T08:09:35.951779Z","shell.execute_reply":"2023-12-17T08:09:39.40429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install --no-index --find-links /tmp/pip/cache/ connected-components-3d","metadata":{"execution":{"iopub.status.busy":"2023-12-17T08:10:17.459573Z","iopub.execute_input":"2023-12-17T08:10:17.460019Z","iopub.status.idle":"2023-12-17T08:10:33.357644Z","shell.execute_reply.started":"2023-12-17T08:10:17.459983Z","shell.execute_reply":"2023-12-17T08:10:33.356606Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Imports","metadata":{}},{"cell_type":"code","source":"import cc3d\n\nimport sys, os\n# for model, ckpt and helper code\nsys.path.append('/kaggle/input/blood-vessel-segmentation-01')\n\nfrom helper import *\n\nimport cv2\nimport pandas as pd\nfrom glob import glob\nimport numpy as np\nfrom skimage.filters import apply_hysteresis_threshold\n\nfrom timeit import default_timer as timer\nimport gc\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nprint('IMPORT OK  !!!!')\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-17T08:12:33.82142Z","iopub.execute_input":"2023-12-17T08:12:33.822138Z","iopub.status.idle":"2023-12-17T08:12:38.990101Z","shell.execute_reply.started":"2023-12-17T08:12:33.822104Z","shell.execute_reply":"2023-12-17T08:12:38.988695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"# for cc3d.dust \ncfg = dotdict( \n    p_threshold=0.29, # best prediction probability threshold \n    cc_connectivity=6, # 6 best connected neighborhood  N.B. only 4,8 (2D) and 26, 18, and 6 (3D) are allowed\n    cc_threshold=6552, # 6552 best  or options from 6474 - 6630 \n    use_tta=True, \n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data ","metadata":{}},{"cell_type":"code","source":"data_dir = \\\n    '/kaggle/input/blood-vessel-segmentation' \nDATA_META = {\n\t'kidney_1_dense': dotdict(\n\t\tname='kidney_1_dense',\n\t\timage_no=(1000, 1000 + 1000),\n\t\timage_dir=f'{data_dir}/train'\n\t),\n\t'kidney_3_dense': dotdict(\n\t\tname='kidney_3_dense',\n\t\timage_no=(496, 996 + 1),\n\t\timage_dir=f'{data_dir}/train'\n\t),\n\t'kidney_5': dotdict(\n\t\tname='kidney_5',\n\t\timage_no=None,\n\t\timage_dir=f'{data_dir}/test'\n\t),\n\t'kidney_6': dotdict(\n\t\tname='kidney_6',\n\t\timage_no=None,\n\t\timage_dir=f'{data_dir}/test'\n\t),\n}\n\n\nmode = 'submit' #'local' \n\nif 'local' in mode:\t\n\tvalid_meta = [DATA_META['kidney_3_dense'], ]\nif 'submit' in mode:\n\tvalid_meta = [DATA_META['kidney_5'], DATA_META['kidney_6']]\n\n\n## io input, etc function\ndef build_file_list(d):\n\tif d.image_no is not None:\n\t\td.file = [f'{d.image_dir}/{d.name.replace(\"kidney_3_dense\",\"kidney_3_sparse\")}/images/{i:04d}.tif' for i in range(*d.image_no)]\n\telse:\n\t\td.file = sorted(glob(f'{d.image_dir}/{d.name}/images/*.tif'))\n\n\nfor d in valid_meta:\n\tbuild_file_list(d)\nprint(valid_meta[0].name, '\\n', valid_meta[0].file[:5]) \nprint('MODE OK  !!!!')","metadata":{"execution":{"iopub.status.busy":"2023-12-17T08:12:53.08882Z","iopub.execute_input":"2023-12-17T08:12:53.08945Z","iopub.status.idle":"2023-12-17T08:12:53.114332Z","shell.execute_reply.started":"2023-12-17T08:12:53.089412Z","shell.execute_reply":"2023-12-17T08:12:53.112645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Functions","metadata":{}},{"cell_type":"code","source":"def file_to_id(f):\n\ts = f.split('/')\n\treturn s[-3] + '_' + s[-1][:-4]\n\n\ndef load_volume(d):\n\tvolume = [\n\t\tcv2.imread(f, cv2.IMREAD_UNCHANGED) for f in d.file\n\t]\n\tvolume = np.stack(volume)\n\treturn volume\n\n\ndef load_truth(d):\n\ttruth = [\n\t\tcv2.imread(f.replace('/images/', '/labels/'), cv2.IMREAD_GRAYSCALE) for f in d.file\n\t]\n\ttruth = np.stack(truth)\n\ttruth = truth // 255\n\treturn truth\n\n\ndef norm_by_percentile(x, low=10, high=99.8, alpha=0.01):\n\txmin = np.percentile(x, low)\n\txmax = np.percentile(x, high)\n\tx = (x - xmin) / (xmax - xmin)\n\tif 1:\n\t\tx[x > 1] = (x[x > 1] - 1) * alpha + 1\n\t\tx[x < 0] = (x[x < 0]) * alpha\n\t# x = np.clip(x,0,1)\n\treturn x\n\n\ndef rle_encode(mask):\n\tpixel = mask.flatten()\n\tpixel = np.concatenate([[0], pixel, [0]])\n\trun = np.where(pixel[1:] != pixel[:-1])[0] + 1\n\trun[1::2] -= run[::2]\n\trle = ' '.join(str(r) for r in run)\n\tif rle == '':\n\t\trle = '1 0'\n\treturn rle\n\n\ndef make_dummy_submission():\n    submission_df = []\n    for d in valid_meta:\n        submission_df.append(\n            pd.DataFrame(data={\n                'id': [file_to_id(f) for f in d.file],\n                'rle': ['1 0'] * len(d.file),\n            })\n        )\n    submission_df = pd.concat(submission_df).reset_index(drop=True)\n    submission_df.to_csv('submission.csv', index=False)\n    return submission_df\n\n\n\nprint('DATASET OK  !!!!')","metadata":{"execution":{"iopub.status.busy":"2023-12-17T08:13:04.818656Z","iopub.execute_input":"2023-12-17T08:13:04.819084Z","iopub.status.idle":"2023-12-17T08:13:04.84268Z","shell.execute_reply.started":"2023-12-17T08:13:04.819051Z","shell.execute_reply":"2023-12-17T08:13:04.841522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model and Submission","metadata":{}},{"cell_type":"code","source":"from model_resnet50d_2d_3d import Net  \n\ncheckpoint_file = \\\n    '/kaggle/input/blood-vessel-segmentation-01/100-resnet50-2d-to-3d-03-00000828.pth'\n \nnet = Net()\nstate_dict = torch.load(checkpoint_file, map_location=lambda storage, loc: storage)['state_dict']\nprint(net.load_state_dict(state_dict, strict=False))  # True\nnet = net.eval()\nnet = net.cuda()\n\n\ndef do_submit():\n    submission_df =[]\n    \n    for d in valid_meta:        \n\n        print('load_volume() ...')\n        volume = load_volume(d)\n        volume = norm_by_percentile(volume)\n        volume = volume.astype(np.float16)\n        D, H, W = volume.shape\n        print(volume.shape)\n\n        # ----\n        depth = 32\n        size  = 640 #640\n        num_D = int(np.ceil((D-8)/(depth-8)))\n        num_H = int(np.ceil((H-32)/(size-32)))\n        num_W = int(np.ceil((W-32)/(size-32)))\n        num_subvolume = num_D*num_H*num_W\n        print('num_D,num_H,num_W', (num_D,num_H,num_W))\n\n        zz = np.linspace(0,D-depth,num_D).astype(int).tolist()\n        yy = np.linspace(0,H-size,num_H).astype(int).tolist()\n        xx = np.linspace(0,W-size,num_W).astype(int).tolist()\n\n\n        prob = np.zeros((D, H, W), dtype=np.float16)\n        prob_count = np.full((D, H, W), fill_value=0.01, dtype=np.float16) \n        start_timer = timer()\n        t=0\n        for z in zz:\n            for y in yy:\n                for x in xx:\n                    print('\\r', f'{t}/{num_subvolume} : {z, y, x} ', time_to_str(timer() - start_timer, 'min'), end='')\n                    t=t+1\n                    image = torch.from_numpy(\n                        volume[z:z + depth, y:y + size, x:x + size]).cuda().unsqueeze(0)\n\n                    vessel = 0\n                    counter = 0\n                    with torch.cuda.amp.autocast(enabled=True):\n                        with torch.no_grad():\n                            v = net(image)\n                            vessel += v\n                            counter += 1\n                            #--- TTA here ----\n                            if cfg.use_tta:\n                                v = net(torch.flip(image,dims=(3,)))\n                                vessel += torch.flip(v,dims=(3,))\n                                counter +=1\n                                v = net(torch.flip(image,dims=(2,)))\n                                vessel += torch.flip(v,dims=(2,))\n                                counter +=1\n                                v = net(torch.flip(image,dims=(1,)))\n                                vessel += torch.flip(v,dims=(1,))\n                                counter +=1\n\n\n                    vessel = vessel/counter\n                    vessel = vessel.half().data.cpu().numpy() #probably memory leak?\n                    vessel = vessel.squeeze(0)\n\n                    prob [z:z + depth, y:y + size, x:x + size] += vessel\n                    prob_count[z:z + depth, y:y + size, x:x + size] += 1\n\n                    #--for debug\n                    if (t<=2) and (mode=='local'):\n                        image = image.float().data.cpu().numpy()\n                        image = image.squeeze(0)\n\n                        m = image.mean(0)\n                        v = vessel.mean(0)\n                        v = np.clip(v*3,0,1)\n                        \n                        plt.figure(figsize=(12,12))\n                        plt.imshow(np.hstack([m, v]),cmap='gray')\n                        plt.show()\n                    #--for debug                    \n                    \n                    del vessel \n                    gc.collect()\n\n        print('')\n        #end of all subvolume \n        prob /= prob_count\n        if (mode == 'local'):\n            np.savez_compressed(f'prob.xyz{d.name}.npz', prob=(prob*255).astype(np.uint8))\n            v = prob.mean(0)\n            v = np.clip(v*10,0,1)**0.5\n            plt.figure(figsize=(12,12))\n            plt.imshow(v,cmap='gray')\n            plt.show()\n          \n        del volume \n        gc.collect()\n        \n        predict = (prob > cfg.p_threshold).astype(np.uint8) \n        # post processing ---\n        if cfg.cc_threshold > 0:\n            predict = cc3d.dust(\n                predict,\n                connectivity=cfg.cc_connectivity,  # 26,\n                threshold=cfg.cc_threshold,\n                in_place=False\n            )\n        \n        submission_df.append(\n            pd.DataFrame(data={\n                'id' : [file_to_id(f) for f in d.file],\n                'rle': [rle_encode(p) for p in predict],\n            })\n        )\n\n        del predict \n        del prob\n        del prob_count\n        gc.collect()\n\n    submission_df = pd.concat(submission_df)\n    submission_df = submission_df.reset_index(drop=True)\n    submission_df.to_csv('submission.csv', index=False)\n    print(submission_df)\n\n\n\nglob_file = glob(f'{data_dir}/test/kidney_5/images/*.tif')\nif (mode == 'submit') and (len(glob_file) == 3):  # cannot do 3d cnn because too few test files\n    submission_df = make_dummy_submission()\n    print(submission_df)\nelse:\n    do_submit()\n\n\nprint('SUBMIT OK!!!')","metadata":{"execution":{"iopub.status.busy":"2023-12-13T06:27:21.013301Z","iopub.execute_input":"2023-12-13T06:27:21.01374Z","iopub.status.idle":"2023-12-13T06:37:20.726622Z","shell.execute_reply.started":"2023-12-13T06:27:21.013701Z","shell.execute_reply":"2023-12-13T06:37:20.725645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Debug for local","metadata":{}},{"cell_type":"code","source":"#debug\nif (mode == 'local'):\n    for d in valid_meta: \n        print(d.name)\n\n        print('load_truth truth ...')\n        truth = load_truth(d)\n        prob  = np.load(f'prob.xyz{d.name}.npz')['prob']/255\n\n        for th in [0.5, 0.4, 0.3, 0.2 ]:\n            predict = (prob > th)\n\n            hit, fp, t_sum, p_sum = np_hit_fp_metric(predict, truth)\n            lb_score = fast_compute_surface_dice_score_from_tensor(predict, truth)\n\n            print(checkpoint_file)\n            print('th=',th)\n            print('hit      :', hit / t_sum)\n            print('fp       :', fp / p_sum)\n            print('lb_score :', lb_score)\n            print('')\n            \n            # post processing ---\n            if cfg.cc_threshold > 0:\n                predict = cc3d.dust(\n                    predict,\n                    connectivity=26,\n                    threshold=cfg.cc_threshold,\n                    in_place=False\n                    )\n                print('after cc3d.dust')\n                hit, fp, t_sum, p_sum = np_hit_fp_metric(predict, truth)\n                lb_score = fast_compute_surface_dice_score_from_tensor(predict, truth)\n                print('hit      :', hit / t_sum)\n                print('fp       :', fp / p_sum)\n                print('lb_score :', lb_score)\n                print('')\n","metadata":{"execution":{"iopub.status.busy":"2023-12-13T06:37:20.728142Z","iopub.execute_input":"2023-12-13T06:37:20.728554Z","iopub.status.idle":"2023-12-13T06:38:23.066255Z","shell.execute_reply.started":"2023-12-13T06:37:20.728517Z","shell.execute_reply":"2023-12-13T06:38:23.0653Z"},"trusted":true},"execution_count":null,"outputs":[]}]}