{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"## Full Pytorch inference code with GPU/CUDA\n## The codes are torch version converted from below public notebooks:\n## https://www.kaggle.com/code/sasrav/inheritance-multiple-network-architectures ### Main code/models - converted to torch.\n## https://www.kaggle.com/code/yoyobar/3d-resnet-baseline-inference ### TTA code.\n## https://www.kaggle.com/code/brettolsen/improving-performance-with-l1-hessian-denoising ### Denoising - converted from cupy to torch.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-14T18:49:17.592639Z","iopub.execute_input":"2023-06-14T18:49:17.593015Z","iopub.status.idle":"2023-06-14T18:49:17.605839Z","shell.execute_reply.started":"2023-06-14T18:49:17.592982Z","shell.execute_reply":"2023-06-14T18:49:17.604874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#install pytorch segmentation models offline\n#https://www.kaggle.com/code/raghaw/offline-install-segmentation-model-pytorch/notebook\n!pip install --no-index --find-links=\"/kaggle/input/segmentation-models-pytorch/\" segmentation-models-pytorch","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:17.645239Z","iopub.execute_input":"2023-06-14T18:49:17.645513Z","iopub.status.idle":"2023-06-14T18:49:34.027975Z","shell.execute_reply.started":"2023-06-14T18:49:17.64549Z","shell.execute_reply":"2023-06-14T18:49:34.02675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import segmentation_models_pytorch as smp","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:34.030368Z","iopub.execute_input":"2023-06-14T18:49:34.030836Z","iopub.status.idle":"2023-06-14T18:49:38.084267Z","shell.execute_reply.started":"2023-06-14T18:49:34.030798Z","shell.execute_reply":"2023-06-14T18:49:38.083285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport math\nimport PIL.Image as Image\nimport torch\nfrom torch.utils.data import DataLoader\nimport torch.nn as nn\nfrom torch.optim import AdamW\nfrom torch.utils.data import DataLoader, Dataset, random_split\nimport torchvision.transforms as T\nimport cv2\nimport albumentations as A\nfrom albumentations.pytorch import ToTensorV2\nimport torch.optim as optim\nfrom torch.optim import lr_scheduler\nimport torch.backends.cudnn as cudnn\nimport torchvision\nimport time\nimport copy\nimport gc\nfrom tqdm.notebook import tqdm\nfrom tqdm import tqdm, trange\nfrom tqdm.auto import tqdm\n\ncudnn.benchmark = True\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nprint(torch.__version__)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:38.085783Z","iopub.execute_input":"2023-06-14T18:49:38.086363Z","iopub.status.idle":"2023-06-14T18:49:39.761661Z","shell.execute_reply.started":"2023-06-14T18:49:38.08633Z","shell.execute_reply":"2023-06-14T18:49:39.760644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#config\nSUBMIT = True #if set to False, evaluate labels accuracy\nTILE_SIZE = 224\nSIZE = TILE_SIZE\nSTRIDE = TILE_SIZE// 2\n\n#specify scan images to load start / end (algorithm will find best positions automatically)\nSTART_SCAN = 28\nEND_SCAN = 36\nINPUT_CHANNELS = END_SCAN - START_SCAN + 1\n\nLR = 0.0001\nLOSS_FUNCTION = torch.nn.BCEWithLogitsLoss()\nTHRESHOLD = 0.45","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.764437Z","iopub.execute_input":"2023-06-14T18:49:39.764793Z","iopub.status.idle":"2023-06-14T18:49:39.7719Z","shell.execute_reply.started":"2023-06-14T18:49:39.764758Z","shell.execute_reply":"2023-06-14T18:49:39.770788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_image(path):\n    images = []\n\n    slideids = range(START_SCAN, END_SCAN + 1)\n    \n    for i in slideids:\n        #read images specified by start_scan and end_scan numbers\n        image = cv2.imread(path + f\"/surface_volume/{i:02}.tif\", 0)\n        \n        pad0 = (TILE_SIZE- image.shape[0] % TILE_SIZE)\n        pad1 = (TILE_SIZE- image.shape[1] % TILE_SIZE)\n        \n        #increase borders to match tile size\n        image = np.pad(image, [(0, pad0), (0, pad1)], constant_values=0)\n        \n        images.append(image)\n        del image\n        gc.collect()\n    #stack images (2D arrays) on axis 2\n    images = np.stack(images, axis=2)\n    \n    return images\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.773639Z","iopub.execute_input":"2023-06-14T18:49:39.77403Z","iopub.status.idle":"2023-06-14T18:49:39.784095Z","shell.execute_reply.started":"2023-06-14T18:49:39.773973Z","shell.execute_reply":"2023-06-14T18:49:39.78316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CustomDataset(Dataset):\n    def __init__(self, images, transform = None):\n        self.images = images\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.images)\n\n    def __getitem__(self, idx):\n        image = self.images[idx]\n        if self.transform:\n            data = self.transform(image = image)\n            image = data['image']\n\n        return image","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.786794Z","iopub.execute_input":"2023-06-14T18:49:39.787565Z","iopub.status.idle":"2023-06-14T18:49:39.799784Z","shell.execute_reply.started":"2023-06-14T18:49:39.787537Z","shell.execute_reply":"2023-06-14T18:49:39.798748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_dataset(paths):\n    \n    images_list = []\n    xyxys = []\n    \n    for path in paths:\n        imageset = read_image(path)\n        print(\"Input images - amount:\" + str(imageset.shape[2]) + \" height: \" + str(imageset.shape[0]) + \" width: \" + str(imageset.shape[1]))\n\n        #split images into multiple smaller images (TILE_SIZE)\n        x1_list = list(range(0, imageset.shape[1] - TILE_SIZE+ 1, TILE_SIZE))\n        y1_list = list(range(0, imageset.shape[0] - TILE_SIZE+ 1, TILE_SIZE))\n\n        for y1 in tqdm(y1_list,desc=\"Building dataset \"):\n            for x1 in x1_list:\n                y2 = y1 + TILE_SIZE\n                x2 = x1 + TILE_SIZE\n                images_list.append(imageset[y1:y2, x1:x2,:])\n                xyxys.append((x1, y1, x2, y2))                 \n                \n    xyxys = np.stack(xyxys)\n    image_dataset = CustomDataset(images = images_list,transform=getTransforms())\n    print(\"Dataset includes \" + str(len(images_list)) + \" - \" + str(TILE_SIZE) + \"x\"  + str(TILE_SIZE) + \" images\")\n\n    return image_dataset, xyxys","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.801646Z","iopub.execute_input":"2023-06-14T18:49:39.802084Z","iopub.status.idle":"2023-06-14T18:49:39.812703Z","shell.execute_reply.started":"2023-06-14T18:49:39.802053Z","shell.execute_reply":"2023-06-14T18:49:39.811824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#trained different pretrained models and selected different architectures to build an average prediction\n#this function is used to create specified models and load the stored weights\ndef createModel(modelnr):\n    if modelnr == 1: #Unet++ with better augmentation\n        model = smp.UnetPlusPlus(\n            encoder_name='resnext50_32x4d', \n            encoder_weights=None, \n            classes=1, \n            in_channels=INPUT_CHANNELS,\n            activation=None)\n        optimizer = optim.Adam(model.parameters(), lr=LR) \n        model = model.to(device)\n        optimizer_ft = optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n        # Decay LR by a factor of 0.1 every 7 EPOCHS\n        exp_lr_scheduler = lr_scheduler.StepLR(optimizer_ft, step_size=7, gamma=0.1)\n        model.load_state_dict(torch.load(\"/kaggle/input/trained-models/best-unetpp-resnext50_32x4d.pth\",map_location=torch.device(device)))\n        \n    elif modelnr == 2: #stride 3 Unet\n        model = smp.Unet(\n            encoder_name='resnext50_32x4d', \n            encoder_weights=None, \n            classes=1, \n            in_channels=INPUT_CHANNELS,\n            activation=None)\n        optimizer = optim.Adam(model.parameters(), lr=LR) \n        model = model.to(device)\n        optimizer_ft = optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n        # Decay LR by a factor of 0.1 every 7 EPOCHS\n        exp_lr_scheduler = lr_scheduler.StepLR(optimizer_ft, step_size=7, gamma=0.1)\n        model.load_state_dict(torch.load(\"/kaggle/input/trained-models/best-unet-resnext50_32x4d.pth\",map_location=torch.device(device)))\n        \n    elif modelnr == 3: #stride 3 Linknet\n        model = smp.Linknet(\n            encoder_name='resnet34', \n            encoder_weights=None, \n            classes=1, \n            in_channels=INPUT_CHANNELS,\n            activation=None)\n        optimizer = optim.Adam(model.parameters(), lr=LR) \n        model = model.to(device)\n        optimizer_ft = optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n        # Decay LR by a factor of 0.1 every 7 EPOCHS\n        exp_lr_scheduler = lr_scheduler.StepLR(optimizer_ft, step_size=7, gamma=0.1)\n        model.load_state_dict(torch.load(\"/kaggle/input/trained-models/best-linknet-resnet34.pth\",map_location=torch.device(device)))\n        \n    elif modelnr == 4: #stride 3 FPN\n        model = smp.FPN(\n            encoder_name='resnet34', \n            encoder_weights=None, \n            classes=1, \n            in_channels=INPUT_CHANNELS,\n            activation=None)\n        optimizer = optim.Adam(model.parameters(), lr=LR) \n        model = model.to(device)\n        optimizer_ft = optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n        # Decay LR by a factor of 0.1 every 7 EPOCHS\n        exp_lr_scheduler = lr_scheduler.StepLR(optimizer_ft, step_size=7, gamma=0.1)\n        model.load_state_dict(torch.load(\"/kaggle/input/trained-models/best-fpn-resnet34.pth\",map_location=torch.device(device)))\n    \n    return model ,optimizer","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.814472Z","iopub.execute_input":"2023-06-14T18:49:39.814724Z","iopub.status.idle":"2023-06-14T18:49:39.831299Z","shell.execute_reply.started":"2023-06-14T18:49:39.814703Z","shell.execute_reply":"2023-06-14T18:49:39.830377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def getTransforms():\n    #to apply normalization\n    augList = [A.Normalize(mean= [0] * INPUT_CHANNELS,std= [1] * INPUT_CHANNELS),ToTensorV2(transpose_mask=True),]\n    aug = A.Compose(augList)\n    return aug","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.832676Z","iopub.execute_input":"2023-06-14T18:49:39.833266Z","iopub.status.idle":"2023-06-14T18:49:39.844041Z","shell.execute_reply.started":"2023-06-14T18:49:39.833235Z","shell.execute_reply":"2023-06-14T18:49:39.843106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def TTA2(x:torch.Tensor,model:nn.Module):\n    #x.shape=(batch,c,h,w)\n    shape=x.shape\n    x=[x,*[torch.rot90(x,k=i,dims=(-2,-1)) for i in range(1,2)]]\n    x=torch.cat(x,dim=0)\n    x=model(x)\n    x=torch.sigmoid(x)\n    x=x.reshape(2,shape[0],*shape[2:])\n    x=[torch.rot90(x[i],k=-i,dims=(-2,-1)) for i in range(2)]\n    x=torch.stack(x,dim=0)\n    return x.mean(0)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.849572Z","iopub.execute_input":"2023-06-14T18:49:39.85099Z","iopub.status.idle":"2023-06-14T18:49:39.85861Z","shell.execute_reply.started":"2023-06-14T18:49:39.850956Z","shell.execute_reply":"2023-06-14T18:49:39.857743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def TTA3(x:torch.Tensor,model:nn.Module):\n    #x.shape=(batch,c,h,w)\n    shape=x.shape\n    x=[x,*[torch.rot90(x,k=i,dims=(-2,-1)) for i in range(1,3)]]\n    x=torch.cat(x,dim=0)\n    x=model(x)\n    x=torch.sigmoid(x)\n    x=x.reshape(3,shape[0],*shape[2:])\n    x=[torch.rot90(x[i],k=-i,dims=(-2,-1)) for i in range(3)]\n    x=torch.stack(x,dim=0)\n    return x.mean(0)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.85975Z","iopub.execute_input":"2023-06-14T18:49:39.861832Z","iopub.status.idle":"2023-06-14T18:49:39.870128Z","shell.execute_reply.started":"2023-06-14T18:49:39.861809Z","shell.execute_reply":"2023-06-14T18:49:39.869116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def TTA4(x:torch.Tensor,model:nn.Module):\n    #x.shape=(batch,c,h,w)\n    shape=x.shape\n    x=[x,*[torch.rot90(x,k=i,dims=(-2,-1)) for i in range(1,4)]]\n    x=torch.cat(x,dim=0)\n    x=model(x)\n    x=torch.sigmoid(x)\n    x=x.reshape(4,shape[0],*shape[2:])\n    x=[torch.rot90(x[i],k=-i,dims=(-2,-1)) for i in range(4)]\n    x=torch.stack(x,dim=0)\n    return x.mean(0)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.871682Z","iopub.execute_input":"2023-06-14T18:49:39.872005Z","iopub.status.idle":"2023-06-14T18:49:39.88169Z","shell.execute_reply.started":"2023-06-14T18:49:39.871976Z","shell.execute_reply":"2023-06-14T18:49:39.880681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def predictTest(model,x):\n    yMask = []\n    #set model to evaluation mode\n    model.eval()\n\n    #turn off gradient tracking\n    with torch.no_grad():\n        for i in tqdm(range(len(x)), desc = \"Predict \" + str(model.name)):\n\n            xPart = x[i]\n            xPart = xPart.unsqueeze(0).cuda()\n            yPred = model(xPart)\n            \n            predMask = yPred.squeeze()\n            predMask = torch.sigmoid(predMask)            \n\n            #append predictions to list\n         #   yMask.append(predMask.cpu().numpy())\n            yMask.append(predMask)\n            \n    #convert list to array\n   # yMask = np.array(yMask)\n   \n    return yMask","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.883109Z","iopub.execute_input":"2023-06-14T18:49:39.883428Z","iopub.status.idle":"2023-06-14T18:49:39.892645Z","shell.execute_reply.started":"2023-06-14T18:49:39.8834Z","shell.execute_reply":"2023-06-14T18:49:39.891681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def predictTestTTA(model,x, no=1):\n    yMask = []\n    #set model to evaluation mode\n    model.eval()\n\n    #turn off gradient tracking\n    with torch.no_grad():\n        for i in tqdm(range(len(x)), desc = \"Predict \" + str(model.name)):\n\n            xPart = x[i]\n            xPart = xPart.unsqueeze(0).cuda()\n          #  yPred = model(xPart)\n            \n          #  predMask = yPred.squeeze()\n        #    predMask = torch.sigmoid(predMask)\n            if no==2:\n                predMask=TTA2(xPart,model)\n            elif no==3:\n                predMask=TTA3(xPart,model)\n            elif no==4:\n                predMask=TTA4(xPart,model)\n\n            #append predictions to list\n           # yMask.append(predMask.cpu().numpy())\n            yMask.append(predMask)\n            \n    #convert list to array\n  #  yMask = np.array(yMask)\n   \n    return yMask","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.894121Z","iopub.execute_input":"2023-06-14T18:49:39.894485Z","iopub.status.idle":"2023-06-14T18:49:39.904397Z","shell.execute_reply.started":"2023-06-14T18:49:39.894456Z","shell.execute_reply":"2023-06-14T18:49:39.903497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Adapted from https://www.kaggle.com/code/stainsby/fast-tested-rle/notebook\n# and https://www.kaggle.com/code/kotaiizuka/faster-rle/notebook\ndef rle(output):\n    #pixels = np.where(output.flatten().cpu() > THRESHOLD, 1, 0).astype(np.uint8)\n    pixels = np.where(output.flatten() > THRESHOLD, 1, 0).astype(np.uint8)\n    pixels[0] = 0\n    pixels[-1] = 0\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 2\n    runs[1::2] = runs[1::2] - runs[:-1:2]\n    return ' '.join(str(x) for x in runs)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.906526Z","iopub.execute_input":"2023-06-14T18:49:39.907012Z","iopub.status.idle":"2023-06-14T18:49:39.917897Z","shell.execute_reply.started":"2023-06-14T18:49:39.906982Z","shell.execute_reply":"2023-06-14T18:49:39.917044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#best accuracy models 6,7,8,9 with different architectures \n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.919489Z","iopub.execute_input":"2023-06-14T18:49:39.919752Z","iopub.status.idle":"2023-06-14T18:49:39.926925Z","shell.execute_reply.started":"2023-06-14T18:49:39.919731Z","shell.execute_reply":"2023-06-14T18:49:39.926097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def preparePlot(predMaskAvg, predMask1, predMask2, predMask3, predMask4, labelMask = None, threshold = 0.45):\n\n    #initialize our figure\n    cols = 6\n    if labelMask is None:\n        cols -= 1\n    figure, ax = plt.subplots(nrows=1, ncols=cols, figsize=(10, 10))\n    \n    \n    #plot \n    ax[0].imshow((predMaskAvg > threshold) * 255)\n    ax[1].imshow((predMask1 > threshold) * 255)\n    ax[2].imshow((predMask2 > threshold) * 255)\n    ax[3].imshow((predMask3 > threshold) * 255)\n    ax[4].imshow((predMask4 > threshold) * 255)\n    if labelMask is not None:\n        ax[5].imshow((labelMask * 255).astype('int8'))\n        \n    #set the titles of the subplots\n    ax[0].set_title(\"Pred Avg \" + str(threshold))\n    ax[1].set_title(\"Unet++ \"  + str(threshold))\n    ax[2].set_title(\"Unet \"  + str(threshold))\n    ax[3].set_title(\"Linknet \"  + str(threshold))\n    ax[4].set_title(\"FPN \"  + str(threshold))\n    \n    if labelMask is not None:\n        ax[5].set_title(\"Label\")\n\n    figure.tight_layout()\n \n    figure.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.928057Z","iopub.execute_input":"2023-06-14T18:49:39.928355Z","iopub.status.idle":"2023-06-14T18:49:39.939683Z","shell.execute_reply.started":"2023-06-14T18:49:39.928333Z","shell.execute_reply":"2023-06-14T18:49:39.939001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def unpackImage(imageList,pad_h,pad_w):\n    #assemble predicted image parts to one big image\n    image = torch.zeros(size=(pad_h, pad_w)).cuda()\n    \n    len_h = pad_h / TILE_SIZE\n    len_w = pad_w / TILE_SIZE\n    \n    h = 0\n    w = 0\n\n    #fill predicted image with partial predictions\n    for i in range(len(imageList)):\n        image[h:h + TILE_SIZE,w:w + TILE_SIZE] = imageList[i,:,:]\n        w += TILE_SIZE\n        if w >= pad_w:\n            h += TILE_SIZE\n            w = 0\n    \n    return image","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.941154Z","iopub.execute_input":"2023-06-14T18:49:39.941797Z","iopub.status.idle":"2023-06-14T18:49:39.953584Z","shell.execute_reply.started":"2023-06-14T18:49:39.941767Z","shell.execute_reply":"2023-06-14T18:49:39.952682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.954978Z","iopub.execute_input":"2023-06-14T18:49:39.955444Z","iopub.status.idle":"2023-06-14T18:49:39.963538Z","shell.execute_reply.started":"2023-06-14T18:49:39.955415Z","shell.execute_reply":"2023-06-14T18:49:39.962643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Attempt to import GPU-accelerated numpy routines\n# try:\n#     import cupy as cp\n#     xp = cp\n# except ImportError:\n#     xp = np\nxp = torch\n\ndelta_lookup = {\n    \"xx\": xp.tensor([[1, -2, 1]], dtype=float).cuda(),\n    \"yy\": xp.tensor([[1], [-2], [1]], dtype=float).cuda(),\n    \"xy\": xp.tensor([[1, -1], [-1, 1]], dtype=float).cuda(),\n}\n\ndef operate_derivative(img_shape, pair):\n    assert len(img_shape) == 2\n    delta = delta_lookup[pair]\n    fft = xp.fft.fftn(delta, img_shape).cuda()\n    return fft * xp.conj(fft).cuda()\n\ndef soft_threshold(vector, threshold):\n    return xp.sign(vector).cuda() * xp.maximum(xp.abs(vector).cuda() - threshold, torch.tensor(0).cuda()).cuda()\n\ndef back_diff(input_image, dim):\n    assert dim in (0, 1)\n  #  r, n = xp.shape(input_image)\n    r, n = input_image.shape\n    size = xp.tensor((r, n)).cuda()\n    position = xp.zeros(2, dtype=int).cuda()\n    temp1 = xp.zeros((r+1, n+1), dtype=float).cuda()\n    temp2 = xp.zeros((r+1, n+1), dtype=float).cuda()\n    \n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] += 1\n    position[dim] += 1\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] -= 1\n    return temp1[0:size[0], 0:size[1]]\n\ndef forward_diff(input_image, dim):\n    assert dim in (0, 1)\n   # r, n = xp.shape(input_image)\n    r, n = input_image.shape\n    size = xp.tensor((r, n)).cuda()\n    position = xp.zeros(2, dtype=int).cuda()\n    temp1 = xp.zeros((r+1, n+1), dtype=float).cuda()\n    temp2 = xp.zeros((r+1, n+1), dtype=float).cuda()\n        \n    size[dim] += 1\n    position[dim] += 1\n\n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] -= 1\n    temp2[0:size[0], 0:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] += 1\n    return -temp1[position[0]:size[0], position[1]:size[1]]\n\ndef iter_deriv(input_image, b, scale, mu, dim1, dim2):\n    g = back_diff(forward_diff(input_image, dim1), dim2)\n    d = soft_threshold(g + b, 1 / mu)\n    b = b + (g - d)\n    L = scale * back_diff(forward_diff(d - b, dim2), dim1)\n    return L, b\n\ndef iter_xx(*args):\n    return iter_deriv(*args, dim1=1, dim2=1)\n\ndef iter_yy(*args):\n    return iter_deriv(*args, dim1=0, dim2=0)\n\ndef iter_xy(*args):\n    return iter_deriv(*args, dim1=0, dim2=1)\n\ndef iter_sparse(input_image, bsparse, scale, mu):\n    d = soft_threshold(input_image + bsparse, 1 / mu)\n    bsparse = bsparse + (input_image - d)\n    Lsparse = scale * (d - bsparse)\n    return Lsparse, bsparse\n\ndef denoise_image(input_image, iter_num=100, fidelity=150, sparsity_scale=10, continuity_scale=0.5, mu=1):\n  #  image_size = xp.shape(input_image)\n    image_size = input_image.shape\n    #print(\"Initialize denoising\")\n    norm_array = (\n        operate_derivative(image_size, \"xx\").cuda() + \n        operate_derivative(image_size, \"yy\").cuda() + \n        2 * operate_derivative(image_size, \"xy\").cuda()\n    )\n    norm_array += (fidelity / mu) + sparsity_scale ** 2\n    b_arrays = {\n        \"xx\": xp.zeros(image_size, dtype=float).cuda(),\n        \"yy\": xp.zeros(image_size, dtype=float).cuda(),\n        \"xy\": xp.zeros(image_size, dtype=float).cuda(),\n        \"L1\": xp.zeros(image_size, dtype=float).cuda(),\n    }\n    g_update = xp.multiply(fidelity / mu, input_image).cuda()\n    for i in tqdm(range(iter_num), total=iter_num):\n        #print(f\"Starting iteration {i+1}\")\n        g_update = xp.fft.fftn(g_update).cuda()\n        if i == 0:\n            g = xp.fft.ifftn(g_update / (fidelity / mu)).real.cuda()\n        else:\n            g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real.cuda()\n        g_update = xp.multiply((fidelity / mu), input_image).cuda()\n        \n        #print(\"XX update\")\n        L, b_arrays[\"xx\"] = iter_xx(g, b_arrays[\"xx\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"YY update\")\n        L, b_arrays[\"yy\"] = iter_yy(g, b_arrays[\"yy\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"XY update\")\n        L, b_arrays[\"xy\"] = iter_xy(g, b_arrays[\"xy\"], 2 * continuity_scale, mu)\n        g_update += L\n        \n        #print(\"L1 update\")\n        L, b_arrays[\"L1\"] = iter_sparse(g, b_arrays[\"L1\"], sparsity_scale, mu)\n        g_update += L\n        \n    g_update = xp.fft.fftn(g_update)\n    g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real.cuda()\n    \n    g[g < 0] = 0\n    g -= g.min()\n    g /= g.max()\n    return g","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:39.965035Z","iopub.execute_input":"2023-06-14T18:49:39.965347Z","iopub.status.idle":"2023-06-14T18:49:42.85854Z","shell.execute_reply.started":"2023-06-14T18:49:39.965318Z","shell.execute_reply":"2023-06-14T18:49:42.857619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def subpredict(path):\n    \n    paths = []\n    paths.append(path)\n    \n    if os.path.isfile(path + \"/inklabels.png\"):\n        labelExist = True\n        labelPath = path + \"/inklabels.png\"\n        label_mask = cv2.imread(labelPath, 0)\n        labelImage = label_mask\n        label_mask = (label_mask / 255).astype('float32')\n    else:\n        label_mask = None\n        label_mask = label_mask\n        labelExist = False\n        labelImage = label_mask\n\n    test_set, _ = create_dataset(paths)\n    \n    model1, _ = createModel(1) #unet++\n    model2, _ = createModel(2) #unet\n    model3, _ = createModel(3) #Linknet\n    model4, _ = createModel(4) #FPN\n\n    yMask1 = predictTestTTA(model1,test_set,no=4) #unet++\n    yMask2 = predictTestTTA(model2,test_set,no=4) #unet\n    yMask3 = predictTestTTA(model3,test_set,no=4) #linknet\n    yMask4 = predictTestTTA(model4,test_set,no=4) #fdn\n    \n    \n    del model1, model2, model3, model4, test_set\n    torch.cuda.empty_cache()\n    gc.collect()\n    \n    return yMask1, yMask2, yMask3, yMask4, labelExist, label_mask, labelImage","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:42.859879Z","iopub.execute_input":"2023-06-14T18:49:42.860305Z","iopub.status.idle":"2023-06-14T18:49:42.869905Z","shell.execute_reply.started":"2023-06-14T18:49:42.860274Z","shell.execute_reply":"2023-06-14T18:49:42.868747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def predictFolderTTA(path):\n    startTime = time.time()\n    #predict for submission\n    TILE_SIZE= 224\n    SIZE= TILE_SIZE\n    \n    #record accuracies\n    acTreshold = []\n    acAvg = []\n    acUnetpp = []\n    acUnet = []\n    acLink = []\n    acFdn = []\n\n    with torch.no_grad():\n        with torch.cuda.amp.autocast(enabled=True):\n            yMask1, yMask2, yMask3, yMask4, labelExist, label_mask, labelImage = subpredict(path)\n            torch.cuda.empty_cache()\n            gc.collect()\n            yMask1 = torch.stack(yMask1)\n            yMask2 = torch.stack(yMask2)\n            yMask3 = torch.stack(yMask3)\n            yMask4 = torch.stack(yMask4)\n            stacked_tensor = torch.stack([yMask1, yMask2, yMask3, yMask4], dim=0)\n            yMask = torch.mean(stacked_tensor, dim=0)\n            del stacked_tensor\n            torch.cuda.empty_cache()\n            gc.collect()\n            #build average prediction  \n           # yMask = (yMask1 + yMask2 + yMask3 + yMask4) /4\n\n            print(\"Predicted mask shape :\" + str(yMask.shape) )\n\n            binary_mask = cv2.imread(path + \"/mask.png\", 0)\n            binary_mask = (binary_mask / 255).astype('float32')\n            binary_mask = torch.tensor(binary_mask).cuda()\n            #original size\n            ori_h = binary_mask.shape[0]\n            ori_w = binary_mask.shape[1]\n\n            #padding again to fit with our tile size\n            pad0 = (TILE_SIZE - binary_mask.shape[0] % TILE_SIZE)\n            pad1 = (TILE_SIZE - binary_mask.shape[1] % TILE_SIZE)\n\n            binary_mask_padded = torch.nn.functional.pad(binary_mask, (0, pad1, 0, pad0), mode='constant', value=0)\n\n            pad_h = binary_mask_padded.shape[0]\n            pad_w = binary_mask_padded.shape[1]\n\n            #assemble predicted image parts to big image\n            yPred = unpackImage(yMask,pad_h,pad_w)\n            yPred1 = unpackImage(yMask1,pad_h,pad_w)\n            yPred2 = unpackImage(yMask2,pad_h,pad_w)\n            yPred3 = unpackImage(yMask3,pad_h,pad_w)\n            yPred4 = unpackImage(yMask4,pad_h,pad_w)\n\n            del yMask1\n            del yMask2\n            del yMask3\n            del yMask4\n\n            gc.collect()\n\n            #back to original size\n            yPred = yPred[:ori_h, :ori_w] \n            yPred1 = yPred1[:ori_h, :ori_w] \n            yPred2 = yPred2[:ori_h, :ori_w]\n            yPred3 = yPred3[:ori_h, :ori_w] \n            yPred4 = yPred4[:ori_h, :ori_w] \n\n            #multiply with mask to set predictions to zero if corresponding mask is zero\n            yPred = torch.multiply(binary_mask, yPred)\n            yPred1 = torch.multiply(binary_mask, yPred1)\n            yPred2 = torch.multiply(binary_mask, yPred2)\n            yPred3 = torch.multiply(binary_mask, yPred3)\n            yPred4 = torch.multiply(binary_mask, yPred4)\n            \n#           #  yPred=xp.array(yPred)\n#             yPred = denoise_image(yPred, iter_num=250)\n#           #  yPred=yPred.get()\n#             yPred = torch.multiply(binary_mask, yPred)\n\n    #show predictions \n    preparePlot(yPred.cpu(), yPred1.cpu(), yPred2.cpu(), yPred3.cpu(), yPred4.cpu(), labelImage, 0.5)\n\n\n    \n    #if labels exist, display accuracy for multiple thresholds\n    \n    if labelExist == True:\n        label_mask = torch.tensor(label_mask).cuda()\n        for t in range(20,70,5):\n            tact = t/100\n            \n            accAvg = (((yPred > tact).float() == label_mask).float().sum() / label_mask.numel()) * 100\n            acc1 = (((yPred1 > tact).float() == label_mask).float().sum() / label_mask.numel()) * 100\n            acc2 = (((yPred2 > tact).float() == label_mask).float().sum() / label_mask.numel()) * 100\n            acc3 = (((yPred3 > tact).float() == label_mask).float().sum() / label_mask.numel()) * 100\n            acc4 = (((yPred4 > tact).float() == label_mask).float().sum() / label_mask.numel()) * 100\n            \n            acTreshold.append(tact)\n            acAvg.append(accAvg.cpu())\n            acUnetpp.append(acc1.cpu())\n            acUnet.append(acc2.cpu())\n            acLink.append(acc3.cpu())\n            acFdn.append(acc4.cpu())           \n            \n            #print(\"Average Accuracy: {:.2f} threshold: {:.2f}\".format(accAvg,tact))\n            #print(\"Accuracy Unet: {:.2f} threshold: {:.2f}\".format(acc2,tact))\n            #print(\"Accuracy Linknet: {:.2f} threshold: {:.2f}\".format(acc3,tact))\n            #print(\"Accuracy FPN: {:.2f} threshold: {:.2f}\".format(acc4,tact))\n            #print(\"*************************************************************\")\n            \n        #plot accuracy \n        plt.figure(2,(15,15))\n        plt.subplot(111)\n        plt.plot(acTreshold, acAvg, label=\"Average\")\n        plt.plot(acTreshold, acUnetpp, label=\"Unet++\")\n        plt.plot(acTreshold, acUnet, label=\"Unet\")\n        plt.plot(acTreshold, acLink, label=\"LinkNet\")\n        plt.plot(acTreshold, acFdn, label=\"FDN\")\n\n        plt.xlabel(\"Activation threshold\")\n        plt.ylabel(\"Accuracy\")\n\n        plt.title(\"Accuracy/Activation\")\n\n        plt.legend()\n\n        plt.show()\n            \n    #multiply with 1 to receive 0,1 values instead of True, False values\n    yPred = (yPred > THRESHOLD) * 1\n    yPred1 = (yPred1 > THRESHOLD) * 1\n    yPred2 = (yPred2 > THRESHOLD) * 1\n    yPred3 = (yPred3 > THRESHOLD) * 1\n    yPred4 = (yPred4 > THRESHOLD) * 1\n\n    #some statistics to compare data structure\n    print(\"yPrediction:\" + \" min: \" + str(yPred.min()) + \" max: \" + str(yPred.max()) + \" shape: \" + str(yPred.shape)  + \" size: \" + str(yPred.size))\n    print(\"Binary Mask:\" + \" min: \" + str(binary_mask.min()) + \" max: \" + str(binary_mask.max()) + \" shape: \" + str(binary_mask.shape) + \" size: \" + str(binary_mask.size))\n    if labelExist == True:\n        print(\"Label Mask:\" + \" min: \" + str(label_mask.min()) + \" max: \" + str(label_mask.max()) + \" shape: \" + str(label_mask.shape) + \" size: \" + str(label_mask.size))\n    print(\"**************************************************************************************************\")\n    #use the average predictions for submission\n    rleOutput = rle(yPred.cpu().numpy())\n\n    endTime = time.time()\n    totalTime = endTime - startTime\n    print(\"Total time: {:.2f} seconds\".format(totalTime))\n    \n\n    \n    return rleOutput","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:42.871468Z","iopub.execute_input":"2023-06-14T18:49:42.871786Z","iopub.status.idle":"2023-06-14T18:49:42.901475Z","shell.execute_reply.started":"2023-06-14T18:49:42.871757Z","shell.execute_reply":"2023-06-14T18:49:42.900627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:42.902935Z","iopub.execute_input":"2023-06-14T18:49:42.903318Z","iopub.status.idle":"2023-06-14T18:49:43.451023Z","shell.execute_reply.started":"2023-06-14T18:49:42.903288Z","shell.execute_reply":"2023-06-14T18:49:43.449194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if SUBMIT == False:\n    _ = predictFolderTTA(\"/kaggle/input/vesuvius-challenge-ink-detection/train/1\")\n    torch.cuda.empty_cache()\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:49.546768Z","iopub.execute_input":"2023-06-14T18:49:49.547131Z","iopub.status.idle":"2023-06-14T18:52:00.839303Z","shell.execute_reply.started":"2023-06-14T18:49:49.547102Z","shell.execute_reply":"2023-06-14T18:52:00.838156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if SUBMIT == False:\n    _ = predictFolderTTA(\"/kaggle/input/vesuvius-challenge-ink-detection/train/2\")\n    torch.cuda.empty_cache()\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:52:37.257839Z","iopub.execute_input":"2023-06-14T18:52:37.258703Z","iopub.status.idle":"2023-06-14T18:57:48.372974Z","shell.execute_reply.started":"2023-06-14T18:52:37.258667Z","shell.execute_reply":"2023-06-14T18:57:48.37196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if SUBMIT == False:\n  #  _ = predictFolderOne(\"/kaggle/input/vesuvius-challenge-ink-detection/train/3\")\n    _ = predictFolderTTA(\"/kaggle/input/vesuvius-challenge-ink-detection/train/3\")\n    torch.cuda.empty_cache()\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:43.457389Z","iopub.status.idle":"2023-06-14T18:49:43.457827Z","shell.execute_reply.started":"2023-06-14T18:49:43.457602Z","shell.execute_reply":"2023-06-14T18:49:43.457623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rleOutput1 = predictFolderTTA(\"/kaggle/input/vesuvius-challenge-ink-detection/test/a\")","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:43.459359Z","iopub.status.idle":"2023-06-14T18:49:43.461487Z","shell.execute_reply.started":"2023-06-14T18:49:43.461293Z","shell.execute_reply":"2023-06-14T18:49:43.46131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rleOutput2 = predictFolderTTA(\"/kaggle/input/vesuvius-challenge-ink-detection/test/b\")","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:43.462806Z","iopub.status.idle":"2023-06-14T18:49:43.46357Z","shell.execute_reply.started":"2023-06-14T18:49:43.463335Z","shell.execute_reply":"2023-06-14T18:49:43.463357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results = []\nresults.append(('a',rleOutput1))\nresults.append(('b',rleOutput2))\nsubmission = pd.DataFrame(results, columns=['Id', 'Predicted'])\nsubmission","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:43.46496Z","iopub.status.idle":"2023-06-14T18:49:43.465678Z","shell.execute_reply.started":"2023-06-14T18:49:43.46545Z","shell.execute_reply":"2023-06-14T18:49:43.465471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2023-06-14T18:49:43.467055Z","iopub.status.idle":"2023-06-14T18:49:43.467792Z","shell.execute_reply.started":"2023-06-14T18:49:43.467559Z","shell.execute_reply":"2023-06-14T18:49:43.46758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}