{"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 torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torchvision\nimport torchvision.transforms as T\nimport torchvision.transforms.functional as TF\nimport cv2\nimport numpy as np\nimport pandas as pd\nfrom torch.utils.data import Dataset, DataLoader\nfrom fastai.vision.all import *\nfrom typing import Optional\nfrom torch.nn.functional import one_hot\nfrom sklearn.model_selection import KFold\nimport random","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:50:47.803456Z","iopub.execute_input":"2024-01-14T13:50:47.804081Z","iopub.status.idle":"2024-01-14T13:50:56.832729Z","shell.execute_reply.started":"2024-01-14T13:50:47.804035Z","shell.execute_reply":"2024-01-14T13:50:56.831245Z"},"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_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_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-14T13:50:56.836249Z","iopub.execute_input":"2024-01-14T13:50:56.836792Z","iopub.status.idle":"2024-01-14T13:50:56.852228Z","shell.execute_reply.started":"2024-01-14T13:50:56.836745Z","shell.execute_reply":"2024-01-14T13:50:56.850987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"volume,target,keys = load_volume('/kaggle/input/blood-vessel-segmentation/train/kidney_1_dense')","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:50:56.854288Z","iopub.execute_input":"2024-01-14T13:50:56.854801Z","iopub.status.idle":"2024-01-14T13:53:14.746253Z","shell.execute_reply.started":"2024-01-14T13:50:56.854759Z","shell.execute_reply":"2024-01-14T13:53:14.744096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n_,Z,X,Y = volume.shape\nZ,X,Y","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:53:14.751536Z","iopub.execute_input":"2024-01-14T13:53:14.752746Z","iopub.status.idle":"2024-01-14T13:53:14.766034Z","shell.execute_reply.started":"2024-01-14T13:53:14.752685Z","shell.execute_reply":"2024-01-14T13:53:14.764775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Given a volume with shape Z*X*Y and the usual XY, ZY and ZX slices you get a total of Z+X+Y images for train. But is simple to triple this data by also considering the slices oblique to those planes. There is a total of two extra sets of slices for each axis, the positive and negative planes. Let's see.","metadata":{}},{"cell_type":"code","source":"class oblique_slices():\n    \n    def __init__(self,volume,axis):\n        SHAPE = list(volume.shape[-3:])\n        SHAPE = SHAPE[:axis]+SHAPE[axis+1:]\n        self.SHAPE = SHAPE\n        self.TOTAL = 2*min(SHAPE) + abs(SHAPE[0]-SHAPE[1]) - 1\n        self.volume = volume\n        self.axis = axis\n        \n    def __len__(self):\n        return self.TOTAL\n        \n    def negative(self,i):\n        H,W = self.SHAPE\n        LINES = min([i+1,min(self.SHAPE),self.TOTAL-i])\n        h = np.arange(i+1)\n        h = h[h < H][-LINES:]\n        w = np.arange(i,-1,-1)\n        w = w[w < W][:LINES]\n        return h,w\n        \n    def positive(self,i):\n        H,W = self.SHAPE\n        LINES = min([i+1,min(self.SHAPE),self.TOTAL-i])\n        h = np.arange(i+1)\n        h = h[h < H][-LINES:]\n        w = np.arange(W-i-1,W)\n        w = w[w > -1][:LINES]\n        return h,w","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:53:14.767581Z","iopub.execute_input":"2024-01-14T13:53:14.767999Z","iopub.status.idle":"2024-01-14T13:53:14.783048Z","shell.execute_reply.started":"2024-01-14T13:53:14.767966Z","shell.execute_reply":"2024-01-14T13:53:14.782106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices = oblique_slices(volume,0)\nlen(slices)","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:53:14.784787Z","iopub.execute_input":"2024-01-14T13:53:14.785215Z","iopub.status.idle":"2024-01-14T13:53:14.805651Z","shell.execute_reply.started":"2024-01-14T13:53:14.785181Z","shell.execute_reply":"2024-01-14T13:53:14.804099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"VOL = np.zeros((Z,X,Y))","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:53:14.807403Z","iopub.execute_input":"2024-01-14T13:53:14.807849Z","iopub.status.idle":"2024-01-14T13:53:14.817125Z","shell.execute_reply.started":"2024-01-14T13:53:14.807786Z","shell.execute_reply":"2024-01-14T13:53:14.815738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import clear_output\n# To visualize how the volume is scanned obliquely, positive planes\nfor axis in [0,1,2]:\n        VOL[:,:,:] = 0\n        slices = oblique_slices(volume,axis)\n        for i in tqdm.tqdm(range(len(slices)//2)):\n            h,w = slices.positive(i)\n            if axis == 0:\n                VOL[:,h,w] = volume[0,:,h,w].T\n            elif axis == 1:\n                VOL[h,:,w] = volume[0,h,:,w]\n            else:\n                VOL[h,w,:] = volume[0,h,w,:]\n                \n            if axis == 0:\n                plt.imshow(VOL[Z//2,:,:])\n            elif axis == 1:\n                plt.imshow(VOL[:,X//2,:])\n            else:\n                plt.imshow(VOL[:,:,Y//2])\n            clear_output(wait=True)\n            plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-14T13:53:14.819224Z","iopub.execute_input":"2024-01-14T13:53:14.821702Z","iopub.status.idle":"2024-01-14T14:27:05.049160Z","shell.execute_reply.started":"2024-01-14T13:53:14.821648Z","shell.execute_reply":"2024-01-14T14:27:05.048139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To visualize how the volume is scanned obliquely, negative planes\nfor axis in [0,1,2]:\n        slices = oblique_slices(volume,axis)\n        VOL[:,:,:] = 0\n        for i in tqdm.tqdm(range(len(slices)//2)):\n            h,w = slices.negative(i)\n            if axis == 0:\n                VOL[:,h,w] = volume[0,:,h,w].T\n            elif axis == 1:\n                VOL[h,:,w] = volume[0,h,:,w]\n            else:\n                VOL[h,w,:] = volume[0,h,w,:]\n            \n            if axis == 0:\n                plt.imshow(VOL[Z//2,:,:])\n            elif axis == 1:\n                plt.imshow(VOL[:,X//2,:])\n            else:\n                plt.imshow(VOL[:,:,Y//2])\n            clear_output(wait=True)\n            plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-14T14:27:05.050778Z","iopub.execute_input":"2024-01-14T14:27:05.051372Z","iopub.status.idle":"2024-01-14T15:00:24.712029Z","shell.execute_reply.started":"2024-01-14T14:27:05.051341Z","shell.execute_reply":"2024-01-14T15:00:24.710708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To visualize the actual slices, positive planes\nfor axis in [0,1,2]:\n        slices = oblique_slices(volume,axis)\n        for i in range(len(slices)//4,len(slices)//2):\n            h,w = slices.positive(i)\n            if axis == 0:\n                plt.imshow(volume[0,:,h,w])\n            elif axis == 1:\n                plt.imshow(volume[0,h,:,w])\n            else:\n                plt.imshow(volume[0,h,w,:])\n            clear_output(wait=True)\n            plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-14T15:00:24.716650Z","iopub.execute_input":"2024-01-14T15:00:24.717161Z","iopub.status.idle":"2024-01-14T15:15:36.210976Z","shell.execute_reply.started":"2024-01-14T15:00:24.717109Z","shell.execute_reply":"2024-01-14T15:15:36.209801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To visualize the actual slices, negative planes\nfor axis in [0,1,2]:\n        slices = oblique_slices(volume,axis)\n        for i in range(len(slices)//4,len(slices)//2):\n            clear_output(wait=True)\n            h,w = slices.negative(i)\n            if axis == 0:\n                plt.imshow(volume[0,:,h,w])\n            elif axis == 1:\n                plt.imshow(volume[0,h,:,w])\n            else:\n                plt.imshow(volume[0,h,w,:])\n            plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-14T15:15:36.212497Z","iopub.execute_input":"2024-01-14T15:15:36.212916Z","iopub.status.idle":"2024-01-14T15:30:26.893601Z","shell.execute_reply.started":"2024-01-14T15:15:36.212885Z","shell.execute_reply":"2024-01-14T15:30:26.892465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"But in oblique slices height spacing is sqrt(2) times larger than width spacing.\nAs augmented data, we could just simply ignore it. But in order to be completely comparable to the rest of the slices we should correct it.","metadata":{}},{"cell_type":"code","source":"def height_correction(img):\n    h,W = img.shape[-2:]\n#   Right now there is h rows space\n#   sqrt(2) times farther than W lines.\n#   We will add the minium amount \n#   of lines to h and revalue them\n#   to fit the W spacing.\n    scale = np.sqrt(2)\n    H = int(scale*h)\n    p = (scale*h - H)/2\n    IMG = np.zeros((H,W))\n    IMG[0,:] = img[0,:]\n    x = p\n    k = 0\n    XL = 0\n    for i in range(1,h):\n#       Just a linear correction between the nearest rows\n#       y = Ax + B\n        A = (img[i,:] - img[i-1,:])/scale\n        B = img[i-1,:]\n        XR = i*scale\n        while x < XR:\n            IMG[k,:]=(A*(x-XL) + B)\n            x += 1\n            k += 1\n            \n        XL = XR\n            \n    return IMG","metadata":{"execution":{"iopub.status.busy":"2024-01-14T15:30:26.895024Z","iopub.execute_input":"2024-01-14T15:30:26.895395Z","iopub.status.idle":"2024-01-14T15:30:26.905817Z","shell.execute_reply.started":"2024-01-14T15:30:26.895365Z","shell.execute_reply":"2024-01-14T15:30:26.904490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"h,w = slices.positive(len(slices)//2)\nimg = volume[0,h,w,:]\nprint('original')\nplt.imshow(img)\nplt.show()\nprint('corrected from original')\nplt.imshow(height_correction(img))\nplt.show()\nimg = (img - np.min(img))/(np.max(img)-np.min(img))\nprint('normalized')\nplt.imshow(img)\nplt.show()\nprint('corrected from normalized')\nplt.imshow(height_correction(img))","metadata":{"execution":{"iopub.status.busy":"2024-01-14T15:30:26.907291Z","iopub.execute_input":"2024-01-14T15:30:26.907615Z","iopub.status.idle":"2024-01-14T15:30:28.616265Z","shell.execute_reply.started":"2024-01-14T15:30:26.907588Z","shell.execute_reply":"2024-01-14T15:30:28.615092Z"},"trusted":true},"execution_count":null,"outputs":[]}]}