{"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":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"That was my first approach. I've trained 2D UNet before, and I found this exercise interesting. However, in my opinion, the memory constraints make this approach less competitive front the 2D one since the model sees less information. I will be glad if it helps someone.","metadata":{}},{"cell_type":"code","source":"import os\nimport math\nimport glob\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","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:48:32.544579Z","iopub.execute_input":"2023-12-19T19:48:32.544918Z","iopub.status.idle":"2023-12-19T19:48:34.534825Z","shell.execute_reply.started":"2023-12-19T19:48:32.544885Z","shell.execute_reply":"2023-12-19T19:48:34.533825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def seed_everything(seed):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = True","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:48:34.536212Z","iopub.execute_input":"2023-12-19T19:48:34.537330Z","iopub.status.idle":"2023-12-19T19:48:34.543873Z","shell.execute_reply.started":"2023-12-19T19:48:34.537293Z","shell.execute_reply":"2023-12-19T19:48:34.543003Z"},"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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-19T19:48:34.545260Z","iopub.execute_input":"2023-12-19T19:48:34.545697Z","iopub.status.idle":"2023-12-19T19:48:34.560575Z","shell.execute_reply.started":"2023-12-19T19:48:34.545602Z","shell.execute_reply":"2023-12-19T19:48:34.559642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/limitz/pytorch-dataset-with-volumetric-augmentations\nclass RandomVolumetricDataset(torch.utils.data.Dataset):\n    ''' Dataset for segmentation of a sparse class. Keeps\n        track of positive samples and favors samples that\n        contain a positive sample.\n        WARNING: do not use in a distributed setting.\n    '''\n    def __init__(self, datasets, PATCH_SIZE=64, length=1000, transform=None):\n        self.volumes = []\n        self.targets = []\n        self.length = length\n        self.shape = [PATCH_SIZE]*3\n        self.transform = transform\n        self.nonzero = []\n        self.means = []\n        self.stds = []\n        self.maxs = []\n        self.mins = []\n        self.target_max = []\n        \n        for dataset in datasets:\n            print(\"loading volume\", dataset)\n            volume, target, _ = load_volume(dataset)\n            self.volumes.append(volume)\n            self.targets.append(target)\n            self.nonzero.append(np.argwhere(target > 0))\n            self.means.append(np.mean(volume))\n #          self.stds.append(np.std(volume))\n            self.maxs.append(np.max(volume))\n            self.mins.append(np.min(volume))\n            self.target_max.append(np.max(target))\n        \n    def __len__(self):\n        return self.length\n\n    def __getitem__(self, idx):\n        vidx = torch.randint(len(self.volumes), (1,)).item()\n        volume = self.volumes[vidx]\n        target = self.targets[vidx]\n        nonzero = self.nonzero[vidx]\n        volume_mean = self.means[vidx]\n#       volume_std = self.stds[vidx]\n        volume_max = self.maxs[vidx]\n        volume_min = self.mins[vidx]\n        target_max = self.target_max[vidx]\n        random = torch.rand(1)\n        \n        if random > 0.9:\n            # Load a random subvolume\n            z,y,x = torch.randint(volume.shape[-3]-self.shape[-3], (1,)).item(), \\\n                    torch.randint(volume.shape[-2]-self.shape[-2], (1,)).item(), \\\n                    torch.randint(volume.shape[-1]-self.shape[-1], (1,)).item()\n        else:\n            # Load a subvolume containing a random sample\n            idx = torch.randint(len(nonzero), (1,)).item()\n            c,z,y,x = nonzero[idx]\n            \n            z += torch.randint(self.shape[-3],(1,)).sub(self.shape[-3]//2).item()\n            y += torch.randint(self.shape[-2],(1,)).sub(self.shape[-2]//2).item()\n            x += torch.randint(self.shape[-1],(1,)).sub(self.shape[-1]//2).item()\n            \n            z = min(max(0,z+self.shape[-3]//2), volume.shape[-3]-self.shape[-3])\n            y = min(max(0,y+self.shape[-2]//2), volume.shape[-2]-self.shape[-2])\n            x = min(max(0,x+self.shape[-3]//2), volume.shape[-1]-self.shape[-1])\n            \n        volume = volume[:,z:z+self.shape[-3], y:y+self.shape[-2], x:x+self.shape[-1]]\n        target = target[:,z:z+self.shape[-3], y:y+self.shape[-2], x:x+self.shape[-1]]\n\n        volume = (volume - volume_min)/(volume_max - volume_min)\n#       volume = volume - volume_mean\n        volume = torch.from_numpy(volume.astype(np.float32))\n        target = torch.from_numpy(target.astype(np.float32))\n        if self.transform is not None:\n            rng = torch.get_rng_state()\n            volume = self.transform(volume)\n            torch.set_rng_state(rng)\n            target = self.transform(target)\n            \n        mask = target > 0#target_max/2\n        target[mask] = 1\n        target[~mask] = 0\n        \n        return volume, target\n                ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:48:34.563831Z","iopub.execute_input":"2023-12-19T19:48:34.564218Z","iopub.status.idle":"2023-12-19T19:48:34.588257Z","shell.execute_reply.started":"2023-12-19T19:48:34.564193Z","shell.execute_reply":"2023-12-19T19:48:34.587218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The augmentations\n# https://www.kaggle.com/code/limitz/pytorch-dataset-with-volumetric-augmentations\nclass  RandomRotationNd(nn.Module):\n    ''' This augmentation first permutes the dimensions as an initial rotation\n        to select the rotation axis, then rotates around the (fixed) z axis. \n        The result is zoomed in to remove empty space and finally permuted \n        once more to move randomize the rotation axis.\n    '''\n    def __init__(self, dims):\n        super().__init__()\n        self.dims = dims\n\n    def forward(self, x):\n        angle = torch.rand(1).item() * 360\n        keep = torch.arange(x.dim() - self.dims)\n        perm = -torch.randperm(self.dims)-1\n        x = x.clone().permute(*[k.item() for k in keep], *[p.item() for p in perm])\n        rad = math.pi * angle / 180\n        scale = abs(math.sin(rad)) + abs(math.cos(rad))\n        for i in range(0, x.shape[-3],8): # presumptuous\n            v = x[...,i:i+8,:,:]\n            w = v.view(-1, *v.shape[-3:])\n            w = TF.rotate(w, angle)#, interpolation=TF.InterpolationMode.BILINEAR)\n            v = w.view(*v.shape)\n            x[...,i:i+8,:,:] = v\n        s = x.shape\n        x = F.interpolate(x, scale_factor=scale, mode=\"bilinear\")\n        x = TF.center_crop(x, s[-2:])\n        perm = -torch.randperm(self.dims)-1\n        x = x.permute(*[k.item() for k in keep], *[p.item() for p in perm])\n        return x\n\nclass RandomRot90Nd(nn.Module):\n    def __init__(self, dims):\n        super().__init__()\n        self.dims = dims\n\n    def forward(self, x):\n        dims = -torch.randperm(self.dims)[:2]-1\n        dims = [d.item() for d in dims]\n        rot = torch.randint(4, (1,)).item()\n        return x.rot90(rot, dims)\n\nclass RandomPermuteNd(nn.Module):\n    def __init__(self, dims):\n        super().__init__()\n        self.dims = dims\n\n    def forward(self, x):\n        perm = -torch.randperm(self.dims)-1\n        keep = torch.arange(x.dim() - self.dims)\n        return x.permute(*[k.item() for k in keep], *[p.item() for p in perm])\n\nclass RandomFlipNd(nn.Module):\n    def  __init__(self, dims, p=0.5):\n        super().__init__()\n        self.dims = dims\n        self.p = p\n        \n    def forward(self, x):\n        for i in range(self.dims):\n            if torch.rand(1) < self.p:\n                x = x.flip(-i-1)\n        return x\n\nclass ToDevice(nn.Module):\n    ''' Sometimes it helps to move the tensor to the gpu before augmentations like\n        rotation. Note however that you need to set num_workers to 0 in the dataloader\n    '''\n    def __init__(self, device=None):\n        super().__init__()\n        self.device = device or (\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n    def forward(self, x):\n        return x.to(self.device)\n ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:48:34.589585Z","iopub.execute_input":"2023-12-19T19:48:34.590721Z","iopub.status.idle":"2023-12-19T19:48:34.611301Z","shell.execute_reply.started":"2023-12-19T19:48:34.590686Z","shell.execute_reply":"2023-12-19T19:48:34.610312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Constructing the dataset\n# https://www.kaggle.com/code/limitz/pytorch-dataset-with-volumetric-augmentations\ntransform = T.Compose((ToDevice(), RandomRotationNd(3), RandomFlipNd(3)))\nds_train = RandomVolumetricDataset([\n    \"/kaggle/input/blood-vessel-segmentation/train/kidney_1_dense\",\n    \"/kaggle/input/blood-vessel-segmentation/train/kidney_3_dense\"\n\n    ], PATCH_SIZE = 64, length=1000, transform=transform)\nds_valid = RandomVolumetricDataset([\n    \"/kaggle/input/blood-vessel-segmentation/train/kidney_2\"\n],PATCH_SIZE = 64, length=500)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:48:34.612492Z","iopub.execute_input":"2023-12-19T19:48:34.612849Z","iopub.status.idle":"2023-12-19T19:54:40.424645Z","shell.execute_reply.started":"2023-12-19T19:48:34.612798Z","shell.execute_reply":"2023-12-19T19:54:40.423754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ds_train.__getitem__(0)[1].shape","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:40.425879Z","iopub.execute_input":"2023-12-19T19:54:40.426206Z","iopub.status.idle":"2023-12-19T19:54:45.293542Z","shell.execute_reply.started":"2023-12-19T19:54:40.426179Z","shell.execute_reply":"2023-12-19T19:54:45.292605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# One problem with this approach.\n# Even making targets more probable to being chosen \n# they are much less present in samples\n# So the model will learn to predict 0\nnp.unique(ds_valid.__getitem__(0)[1],return_counts=True)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.295078Z","iopub.execute_input":"2023-12-19T19:54:45.295439Z","iopub.status.idle":"2023-12-19T19:54:45.313713Z","shell.execute_reply.started":"2023-12-19T19:54:45.295405Z","shell.execute_reply":"2023-12-19T19:54:45.312821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\ndef double_convolution(in_channels, mid_channels, out_channels):\n    \"\"\"\n    In the original paper implementation, the convolution operations were\n    not padded but we are padding them here. This is because, we need the \n    output result size to be same as input size.\n    \"\"\"\n    conv_op = nn.Sequential(\n                  nn.Conv3d(in_channels, mid_channels, kernel_size=3, padding=1, bias=False),\n                  nn.BatchNorm3d(mid_channels, eps=1e-4),\n                  nn.ReLU(inplace=True),\n                  nn.Conv3d(mid_channels, out_channels, kernel_size=3, padding=1, bias=False),\n                  nn.BatchNorm3d(out_channels, eps=1e-4),\n                  nn.ReLU(inplace=True)\n                            )\n    \n    return conv_op","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.315144Z","iopub.execute_input":"2023-12-19T19:54:45.315693Z","iopub.status.idle":"2023-12-19T19:54:45.323251Z","shell.execute_reply.started":"2023-12-19T19:54:45.315665Z","shell.execute_reply":"2023-12-19T19:54:45.322355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"device = 'cuda' if torch.cuda.is_available() else 'cpu'\n\nclass UNet3D(nn.Module):\n#   I set wrapped layers to cuda at class definition because that implementation,\n#   although I like it for clean, don't allow changes on wrapped layers\n#   It doesn't allow model.state_dict() neither, you must save all the model\n    def __init__(self, num_classes):\n        super(UNet3D, self).__init__()\n        LAYERS = [32, 64, 128, 256, 320, 320]\n        \n        self.LEFT = [double_convolution(1, LAYERS[0], LAYERS[0]).to(device)]\n        for i in range(1,len(LAYERS)):\n            self.LEFT += [double_convolution(LAYERS[i-1], LAYERS[i], LAYERS[i]).to(device)]\n            \n        self.RIGHT = []\n        for i in range(len(LAYERS)-1):\n            self.RIGHT += [double_convolution(LAYERS[i]*2, LAYERS[i]*2, LAYERS[i]).to(device)]\n            \n        self.DOWN = []\n#       Those may be should better be MaxPool3d\n        for i in range(len(LAYERS)-1):\n            self.DOWN += [\n                nn.Conv3d(\n                    in_channels=LAYERS[i],\n                    out_channels=LAYERS[i],\n                    kernel_size=2,\n                    stride=2).to(device)]\n            \n        self.UP = []\n        for i in range(len(LAYERS)-1):\n            self.UP += [\n                nn.ConvTranspose3d(\n                    in_channels=LAYERS[i+1],\n                    out_channels=LAYERS[i],\n                    kernel_size=2,\n                    stride=2).to(device)]\n            \n        # output => `out_channels` as per the number of classes.\n        self.OUT = nn.Conv3d(\n            in_channels=LAYERS[0], out_channels=num_classes, \n            kernel_size=1\n        ).to(device) \n            \n    def forward(self, x):\n        S = []\n        for i in range(len(self.LEFT)-1):\n            x = self.LEFT[i](x)\n            S += [x]\n            x = self.DOWN[i](x)\n            \n        # BRIDGE\n        x = self.LEFT[-1](x)\n        \n        for i in range(len(self.LEFT)-2,-1,-1):\n            x = self.RIGHT[i](torch.cat((S[i],self.UP[i](x)),1))\n            \n        x = self.OUT(x)\n            \n        return x","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.324579Z","iopub.execute_input":"2023-12-19T19:54:45.324910Z","iopub.status.idle":"2023-12-19T19:54:45.341438Z","shell.execute_reply.started":"2023-12-19T19:54:45.324885Z","shell.execute_reply":"2023-12-19T19:54:45.340593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from typing import Optional\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom torch.nn.functional import one_hot\n\n\n# based on:\n# https://github.com/kevinzakka/pytorch-goodies/blob/master/losses.py\n\nclass myDiceLoss(nn.Module):\n    r\"\"\"Criterion that computes Sørensen-Dice Coefficient loss.\n\n    According to [1], we compute the Sørensen-Dice Coefficient as follows:\n\n    .. math::\n\n        \\text{Dice}(x, class) = \\frac{2 |X| \\cap |Y|}{|X| + |Y|}\n\n    where:\n       - :math:`X` expects to be the scores of each class.\n       - :math:`Y` expects to be the one-hot tensor with the class labels.\n\n    the loss, is finally computed as:\n\n    .. math::\n\n        \\text{loss}(x, class) = 1 - \\text{Dice}(x, class)\n\n    [1] https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient\n\n    Shape:\n        - Input: :math:`(N, C, H, W)` where C = number of classes.\n        - Target: :math:`(N, H, W)` where each value is\n          :math:`0 ≤ targets[i] ≤ C−1`.\n\n    Examples:\n        >>> N = 5  # num_classes\n        >>> loss = tgm.losses.DiceLoss()\n        >>> input = torch.randn(1, N, 3, 5, requires_grad=True)\n        >>> target = torch.empty(1, 3, 5, dtype=torch.long).random_(N)\n        >>> output = loss(input, target)\n        >>> output.backward()\n    \"\"\"\n\n    def __init__(self) -> None:\n        super(myDiceLoss, self).__init__()\n        self.eps: float = 1e-6\n\n    def forward(\n            self,\n            input: torch.Tensor,\n            target: torch.Tensor) -> torch.Tensor:\n        if not torch.is_tensor(input):\n            raise TypeError(\"Input type is not a torch.Tensor. Got {}\"\n                            .format(type(input)))\n        if not len(input.shape) == 5:\n            raise ValueError(\"Invalid input shape, we expect BxNxDxHxW. Got: {}\"\n                             .format(input.shape))\n        if not input.shape[-3:] == target.shape[-3:]:\n            raise ValueError(\"input and target shapes must be the same. Got: {}\"\n                             .format(input.shape, input.shape))\n        if not input.device == target.device:\n            raise ValueError(\n                \"input and target must be in the same device. Got: {}\" .format(\n                    input.device, target.device))\n        # compute softmax over the classes axis\n        input_soft = F.softmax(input, dim=1)\n\n        # create the labels one hot tensor\n        target_one_hot = F.one_hot(target.long(), num_classes=input.shape[1])\n        target_one_hot = torch.swapaxes(target_one_hot,1,-1).squeeze(-1)\n\n        # compute the actual dice score\n        dims = (1, 2, 3, 4)\n        intersection = torch.sum(input_soft * target_one_hot, dims)\n        cardinality = torch.sum(input_soft + target_one_hot, dims)\n\n        dice_score = 2. * intersection / (cardinality + self.eps)\n        return torch.mean(1. - dice_score)\n\n\n\n######################\n# functional interface\n######################\n\n\ndef dice_loss(\n        input: torch.Tensor,\n        target: torch.Tensor) -> torch.Tensor:\n    r\"\"\"Function that computes Sørensen-Dice Coefficient loss.\n\n    See :class:`~torchgeometry.losses.DiceLoss` for details.\n    \"\"\"\n    return myDiceLoss()(input, target)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.342705Z","iopub.execute_input":"2023-12-19T19:54:45.343005Z","iopub.status.idle":"2023-12-19T19:54:45.357441Z","shell.execute_reply.started":"2023-12-19T19:54:45.342972Z","shell.execute_reply":"2023-12-19T19:54:45.356504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_batch_size = 4\nvalid_batch_size = 4\ndl_train = torch.utils.data.DataLoader(ds_train,batch_size=train_batch_size,num_workers=0)\ndl_valid = torch.utils.data.DataLoader(ds_valid,batch_size=valid_batch_size,num_workers=0)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.358590Z","iopub.execute_input":"2023-12-19T19:54:45.358878Z","iopub.status.idle":"2023-12-19T19:54:45.369775Z","shell.execute_reply.started":"2023-12-19T19:54:45.358854Z","shell.execute_reply":"2023-12-19T19:54:45.369020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport os, gc\nimport numpy as np\nfrom sklearn.model_selection import KFold\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom fastai.vision.all import *\nseed_everything(2023)\nmodel = UNet3D(2)\ndls = DataLoaders(dl_train,dl_valid)\nlearn = Learner(dls, model, loss_func=dice_loss,cbs=[GradientClip(3.0)])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:45.373541Z","iopub.execute_input":"2023-12-19T19:54:45.373815Z","iopub.status.idle":"2023-12-19T19:54:50.259191Z","shell.execute_reply.started":"2023-12-19T19:54:45.373790Z","shell.execute_reply":"2023-12-19T19:54:50.258352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"learn.fit_one_cycle(30)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:54:50.260353Z","iopub.execute_input":"2023-12-19T19:54:50.260716Z","iopub.status.idle":"2023-12-19T21:21:53.185166Z","shell.execute_reply.started":"2023-12-19T19:54:50.260687Z","shell.execute_reply":"2023-12-19T21:21:53.184231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"torch.save(model,'UNet3D_64')\ndel model","metadata":{"execution":{"iopub.status.busy":"2023-12-19T21:21:54.728746Z","iopub.execute_input":"2023-12-19T21:21:54.729544Z","iopub.status.idle":"2023-12-19T21:21:55.029176Z","shell.execute_reply.started":"2023-12-19T21:21:54.729508Z","shell.execute_reply":"2023-12-19T21:21:55.028148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Inference, just because I've been asked to. But remember that right now the model achieves a validation loss of around 0.15 with 0 predictions only because the proportion of 0 labels in the samples is very imbalanced. We can expect similar volumes for test data sets. Which makes it impossible to process the entire volume at once. So the simplest algorithm we can use will be a patch scan around everywhere.","metadata":{}},{"cell_type":"code","source":"# https://www.kaggle.com/code/limitz/pytorch-dataset-with-volumetric-augmentations\ndef load_pad_volume(dataset, labeled=True, slice_range=None, P=64):\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    D = 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            H,W = image.shape[-2:]\n            volume = np.zeros((1,P*(D//P + 1),\n                                 P*(H//P + 1),\n                                 P*(W//P + 1)), dtype=np.uint16)\n        volume[:,z-offset,:H,:W] = image\n        \n    SHAPE = np.array([D,H,W])\n    \n    return volume, target, keys, SHAPE\n\n# ref.: https://www.kaggle.com/stainsby/fast-tested-rle\ndef rle_encode(img):\n    '''\n    img: numpy array, 1 - mask, 0 - background\n    Returns run length as string formated\n    '''\n    pixels = img.flatten()\n    pixels = np.concatenate([[0], pixels, [0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1\n    runs[1::2] -= runs[::2]\n    \n    mystr = ' '.join(str(x) for x in runs)\n    if mystr == \"\":\n        mystr = \"1 0\"\n    return mystr","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = torch.load('UNet3D_64').to(device)\nP = 64","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Datasets = glob.glob(\"/kaggle/input/blood-vessel-segmentation/test/*\")\nDatasets","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def scan3D(Dataset,P=64):\n    volume, _, keys, SHAPE = load_pad_volume(dataset, labeled=False, slice_range=None)\n    PREDS = torch.zeros(volume.shape[-3:],dtype=torch.int).to(device)\n    STEPS = SHAPE//P + 1\n    \n    for i in range(STEPS[-3]):\n        D = i*P\n        for j in range(STEPS[-2]):\n            H = j*P\n            for j in range(STEPS[-1]):\n                W = k*P\n                SV = volume[...,D:D+P,H:H+P,W:W+P].unsqueeze(0).to(device)\n                PREDS[D:D+P,H:H+P,W:W+P] = torch.argmax(model(SV),axis=-4)\n    \n    D,H,W = SHAPE\n\n    return keys, PREDS[:D,:H,:W]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"All_keys = []\nrle = []\nwith torch.no_grad():\n    for dataset in Datasets:\n        print(dataset)\n        keys, PREDS = load_volume(dataset, labeled=False, slice_range=None)\n        All_keys += keys\n        for PRED in PREDS:\n            rle += [rle_encode(PRED.cpu())]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame({'id':All_keys, 'rle':rle})\nsubmission.head()","metadata":{},"execution_count":null,"outputs":[]}]}