{"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":"none","dataSources":[{"sourceId":61446,"databundleVersionId":6962461,"sourceType":"competition"}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport math\nimport glob\nimport gc\nimport tqdm\nimport cv2\nimport numpy as np\nimport pandas as pd","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-01-12T13:14:56.333371Z","iopub.execute_input":"2024-01-12T13:14:56.333968Z","iopub.status.idle":"2024-01-12T13:14:56.341130Z","shell.execute_reply.started":"2024-01-12T13:14:56.333920Z","shell.execute_reply":"2024-01-12T13:14:56.339839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/limitz/pytorch-dataset-with-volumetric-augmentations\ndef load_volume(dataset, labeled=True, slice_range=None):\n    ''' Load slices into a volume. Keeps the memory requirement\n        as low as possible by using uint8 and uint16 in CPU memory.\n    '''\n    if labeled:\n        path = os.path.join(dataset, \"labels\", \"*.tif\")\n    else:\n        path = os.path.join(dataset, \"images\", \"*.tif\")\n        \n    dataset = sorted(glob.glob(path))\n    volume = None\n    target = None\n    keys = []\n    offset = 0 if slice_range is None else slice_range[0]\n    depth = len(dataset) if slice_range is None else slice_range[1]-slice_range[0]\n    \n    for z, path in enumerate(tqdm.tqdm(dataset)):\n        if slice_range is not None:\n            if z < slice_range[0]: continue\n            if z >= slice_range[1]: continue\n        \n        parts = path.split(os.path.sep)\n        key = parts[-3] + \"_\" + parts[-1].split(\".\")[0]\n        keys.append(key)\n                \n        if labeled:\n            label = cv2.imread(path, cv2.IMREAD_GRAYSCALE)#IMREAD_ANYDEPTH)\n            label = np.array(label,dtype=np.uint8)\n            if target is None:\n                target = np.zeros((1,depth, *label.shape[-2:]), dtype=np.uint8)\n            target[:,z-offset] = label\n        \n        path = path.replace(\"labels\",\"images\")\n        path = path.replace(\"kidney_3_dense\",\"kidney_3_sparse\")\n        image = cv2.imread(path, cv2.IMREAD_GRAYSCALE)#IMREAD_ANYDEPTH)\n        image = np.array(image,dtype=np.uint16)\n        \n        if volume is None:\n            volume = np.zeros((1,depth, *image.shape[-2:]), dtype=np.uint16)\n        volume[:,z-offset] = image\n    \n    return volume, target, keys","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:14:58.811770Z","iopub.execute_input":"2024-01-12T13:14:58.812182Z","iopub.status.idle":"2024-01-12T13:14:58.826792Z","shell.execute_reply.started":"2024-01-12T13:14:58.812145Z","shell.execute_reply":"2024-01-12T13:14:58.825501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"volume,target,_ = load_volume('/kaggle/input/blood-vessel-segmentation/train/kidney_1_dense')","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:15:02.859103Z","iopub.execute_input":"2024-01-12T13:15:02.859480Z","iopub.status.idle":"2024-01-12T13:17:33.811344Z","shell.execute_reply.started":"2024-01-12T13:15:02.859450Z","shell.execute_reply":"2024-01-12T13:17:33.810122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import clear_output\nimport matplotlib.pyplot as plt\nfor vol in volume[0]:\n    pmin,pmax = np.percentile(vol, (1, 99))\n    f, axarr = plt.subplots(1,3)\n    values,counts = np.unique(vol,return_counts=True)\n    axarr[2].plot(values,counts,'b.')\n    axarr[0].imshow(vol)\n    clip_vol = vol.copy()\n    clip_vol[vol < pmin] = pmin\n    clip_vol[vol > pmax] = pmax\n    axarr[1].imshow((clip_vol-pmin)/(pmax - pmin))\n    clear_output(wait=True)\n    plt.show()\n    break","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:18:05.323752Z","iopub.execute_input":"2024-01-12T13:18:05.324946Z","iopub.status.idle":"2024-01-12T13:18:06.185135Z","shell.execute_reply.started":"2024-01-12T13:18:05.324882Z","shell.execute_reply":"2024-01-12T13:18:06.183618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    f, axarr = plt.subplots(1,3)\n    vol = volume[0,-1]\n    values,counts = np.unique(vol,return_counts=True)\n    axarr[2].plot(values,counts,'b.')\n    axarr[0].imshow(vol)\n    clip_vol = vol.copy()\n    clip_vol[vol < pmin] = pmin\n    clip_vol[vol > pmax] = pmax\n    axarr[1].imshow((clip_vol-pmin)/(pmax - pmin))\n    clear_output(wait=True)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:18:11.379756Z","iopub.execute_input":"2024-01-12T13:18:11.380316Z","iopub.status.idle":"2024-01-12T13:18:12.203471Z","shell.execute_reply.started":"2024-01-12T13:18:11.380271Z","shell.execute_reply":"2024-01-12T13:18:12.202148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(\n    len(volume[0])):\n    values,counts = np.unique(volume[0,i,:,:],return_counts=True)\n    clear_output(wait=True)\n    plt.plot(values,counts,'b.')\n    plt.show()\n    break","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:18:15.573757Z","iopub.execute_input":"2024-01-12T13:18:15.574298Z","iopub.status.idle":"2024-01-12T13:18:15.870208Z","shell.execute_reply.started":"2024-01-12T13:18:15.574256Z","shell.execute_reply":"2024-01-12T13:18:15.868830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks like the first peack is basically background","metadata":{}},{"cell_type":"code","source":"median = np.percentile(volume,50)\nmedian","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:18:18.875305Z","iopub.execute_input":"2024-01-12T13:18:18.875850Z","iopub.status.idle":"2024-01-12T13:18:39.295988Z","shell.execute_reply.started":"2024-01-12T13:18:18.875807Z","shell.execute_reply":"2024-01-12T13:18:39.294301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://stackoverflow.com/questions/19206332/gaussian-fit-for-python\nimport pylab as plb\nimport matplotlib.pyplot as plt\nfrom scipy.optimize import curve_fit\nfrom scipy import asarray as ar,exp\n\nvalues,counts = np.unique(volume[0],return_counts=True)\ny = counts[values > median]\nx = values[values > median]\ni = 0\nwhile(y[i+1]<y[i]):i+=1\ny = y[i:]\nx = x[i:]\n\nmean = sum(x*y)/np.sum(y)\nsigma = sum(y*(x-mean)**2)/np.sum(y)\n\ndef gaus(x,a,x0,sigma):\n    return a*np.exp(-(x-x0)**2/(2*sigma**2))\n\npopt_kidney,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma],maxfev=999999)\n\ny = counts\nx = values\n\nplt.plot(x,y,'b+:',label='data')\nplt.plot(x,gaus(x,*popt_kidney),'ro:',label='fit')\nplt.legend()\nplt.title('kidney peack')\nplt.xlabel('values')\nplt.ylabel('counts')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:18:39.298933Z","iopub.execute_input":"2024-01-12T13:18:39.299382Z","iopub.status.idle":"2024-01-12T13:20:12.901822Z","shell.execute_reply.started":"2024-01-12T13:18:39.299346Z","shell.execute_reply":"2024-01-12T13:20:12.900380Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = counts[values < median]\nx = values[values < median]\n\nmean = sum(x*y)/np.sum(y)\nsigma = sum(y*(x-mean)**2)/np.sum(y)\n\npopt_background,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma],maxfev=999999)\n\ny = counts\nx = values\n\nplt.plot(x,y,'b+:',label='data')\nplt.plot(x,gaus(x,*popt_background),'ro:',label='fit')\nplt.legend()\nplt.title('background peack')\nplt.xlabel('values')\nplt.ylabel('counts')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:20:12.903736Z","iopub.execute_input":"2024-01-12T13:20:12.904265Z","iopub.status.idle":"2024-01-12T13:20:13.275523Z","shell.execute_reply.started":"2024-01-12T13:20:12.904222Z","shell.execute_reply":"2024-01-12T13:20:13.274402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We've fitted a gaussian function to each peack\nprint('background peack:',popt_background)\nprint('kidney peack:',popt_kidney)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:20:13.278566Z","iopub.execute_input":"2024-01-12T13:20:13.279920Z","iopub.status.idle":"2024-01-12T13:20:13.287931Z","shell.execute_reply.started":"2024-01-12T13:20:13.279872Z","shell.execute_reply":"2024-01-12T13:20:13.286531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://stackoverflow.com/questions/28766692/how-to-find-the-intersection-of-two-graphs\nx = np.arange(255)\nf = gaus(np.arange(255),*popt_background)\ng = gaus(np.arange(255),*popt_kidney)\n\nplt.plot(x, f, '-')\nplt.plot(x, g, '-')\n\nidx = np.argwhere(np.diff(np.sign(f - g))).flatten()\nplt.plot(x[idx], f[idx], 'ro')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:20:13.289947Z","iopub.execute_input":"2024-01-12T13:20:13.290465Z","iopub.status.idle":"2024-01-12T13:20:13.595737Z","shell.execute_reply.started":"2024-01-12T13:20:13.290423Z","shell.execute_reply":"2024-01-12T13:20:13.594413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idx","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:20:13.597227Z","iopub.execute_input":"2024-01-12T13:20:13.597635Z","iopub.status.idle":"2024-01-12T13:20:13.608125Z","shell.execute_reply.started":"2024-01-12T13:20:13.597599Z","shell.execute_reply":"2024-01-12T13:20:13.606523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"background = volume[0] < idx[1]\nfor axis in [0,1,2]:\n    plt.imshow(background.sum(axis))\n    plt.show()\ndel background","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(len(volume[0])//2):\n    f, axarr = plt.subplots(1,2)\n    img = volume[0,i].copy()\n    img[img < idx[1]] = idx[1]\n    img = (img - idx[1])/(255 - idx[1])\n    axarr[0].imshow(img)\n    axarr[1].imshow(target[0,i])\n    clear_output(wait=True)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:20:13.609945Z","iopub.execute_input":"2024-01-12T13:20:13.610482Z","iopub.status.idle":"2024-01-12T13:42:03.859485Z","shell.execute_reply.started":"2024-01-12T13:20:13.610436Z","shell.execute_reply":"2024-01-12T13:42:03.858134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It works surprisingly well, however it also masks indiscriminately any interior cavity.\nWe can make it better.","metadata":{}},{"cell_type":"code","source":"# https://pypi.org/project/connected-components-3d/\n!pip install connected-components-3d\nimport cc3d\n\n# Get a labeling of the k largest objects in the image.\n# The output will be relabeled from 1 to N.\nbackground= cc3d.largest_k(\n  volume[0] < idx[1], k=1, \n  connectivity=26, delta=0,\n  return_N=False,\n)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:42:03.861450Z","iopub.execute_input":"2024-01-12T13:42:03.861832Z","iopub.status.idle":"2024-01-12T13:44:10.087963Z","shell.execute_reply.started":"2024-01-12T13:42:03.861794Z","shell.execute_reply":"2024-01-12T13:44:10.086570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for axis in [0,1,2]:\n    plt.imshow(background.sum(axis))\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T13:44:10.089944Z","iopub.execute_input":"2024-01-12T13:44:10.090605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(len(volume[0])//2):\n    f, axarr = plt.subplots(1,3)\n    img = volume[0,i].copy()\n    img[background[i]>0] = idx[1]\n    axarr[0].imshow(~(background[i]>0))\n    axarr[1].imshow(img)\n    axarr[2].imshow(target[0,i])\n    clear_output(wait=True)\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DISCLAIMER: At the moment it only works with kidney_1_dense and kidney_3 (if I remember correctly) but fails with kidney_1_voi and kidney_2. Because there is not enough differentiable background peack. May be you're tempted to use it with training, but nothing grants that this can be applied to the test.\nAnd for my experiments is not necessary to remove the background. I did this as excercise and with other possible use than trainig in mind.","metadata":{}}]}