{"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":"import gc \nimport torch \nimport numpy as np \nimport pandas as pd \nfrom pathlib import Path \nimport PIL.Image as Image \nimport matplotlib.pyplot as plt\n\nimport glob \nimport torch.nn as nn \nfrom tqdm import tqdm \nimport torch.utils.data as data\n\nBUFFER_1 = 30 \nWINDOW_1 = 1 \nBUFFER_2 = 60 \nWINDOW_2 = 2 \nBUFFER_3 = 120\nWINDOW_3 = 3\nBATCH_SIZE = 32 \nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nfile_path = '/kaggle/input/vesuvius-challenge-ink-detection/test'","metadata":{"execution":{"iopub.status.busy":"2023-06-14T15:16:30.244871Z","iopub.execute_input":"2023-06-14T15:16:30.245265Z","iopub.status.idle":"2023-06-14T15:16:33.806213Z","shell.execute_reply.started":"2023-06-14T15:16:30.245230Z","shell.execute_reply":"2023-06-14T15:16:33.805000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SubvolumeDataset(data.Dataset):\n    def __init__(self, image_stack, pixels, buffer):\n        self.image_stack = image_stack\n        self.pixels = pixels\n        self.buffer = buffer\n    def __len__(self):\n        return len(self.pixels)\n    def __getitem__(self, index):\n        y, x = self.pixels[index]\n        subvolume = self.image_stack[:, y-self.buffer:y+self.buffer+1, x-self.buffer:x+self.buffer+1].view(1, len(self.image_stack), self.buffer*2+1, self.buffer*2+1)\n        return subvolume\n    \n# 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(img):\n    pixels = img.flatten()\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)\n\n# from io import StringIO\n# def rle_2(img): \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#     f = StringIO() \n#     np.savetxt(f, runs.reshape(1, -1), delimiter=\" \", fmt=\"%d\") \n#     predicted = f.getvalue().strip() \n#     return predicted","metadata":{"_uuid":"7cc9fc6e-c06c-4041-8783-fd50e1fac413","_cell_guid":"8d145586-adc4-43d7-9b51-9188708c127f","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-14T15:16:33.809992Z","iopub.execute_input":"2023-06-14T15:16:33.810642Z","iopub.status.idle":"2023-06-14T15:16:33.828972Z","shell.execute_reply.started":"2023-06-14T15:16:33.810596Z","shell.execute_reply":"2023-06-14T15:16:33.827445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def slice_finder(PREFIX):\n    mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n    mask[:int(mask.shape[0]*.35),:] = 0\n    mask[int(mask.shape[0]*.65):,:] = 0\n    mask[:,:int(mask.shape[1]*.35)] = 0\n    mask[:,int(mask.shape[1]*.65):] = 0\n\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER_1:mask.shape[0]-BUFFER_1, BUFFER_1:mask.shape[1]-BUFFER_1] = True\n\n    grid = np.zeros(mask.shape)\n    i = WINDOW_1\n    while i < grid.shape[0]:\n        j = WINDOW_1\n        while j < grid.shape[1]:\n            grid[i,j] = 1\n            j += 2*WINDOW_1+1\n        i += 2*WINDOW_1+1\n\n    arr_mask = np.array(mask) * not_border * grid    \n    full_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    all_pixels = np.argwhere(full_rect)\n\n    del not_border, arr_mask, grid, full_rect\n    gc.collect()\n\n    highest_avg = 0\n    slices = []\n    \n    k=10\n    while k < 31:\n        image_selections = [k, k+2, k+4, k+6, k+8, k+10, k+12, k+14, k+16, k+18]\n\n        Z_DIM = len(image_selections) # We add x to this if we're including prediction layer(s).\n        img_layers = sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))\n\n        images = [np.array(Image.open(img_layers[i-1]), dtype=np.float32)/65535.0 for i in image_selections]\n        image_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n        del images\n\n        eval_dataset = SubvolumeDataset(image_stack, all_pixels, BUFFER_1)\n\n        del image_stack\n        gc.collect()\n\n        eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\n        del eval_dataset\n\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L1.pth', map_location=DEVICE).to(DEVICE)\n\n        output=[]\n        model.eval()\n        with torch.no_grad():\n            for i, (subvolumes) in enumerate(eval_loader):\n                for j, value in enumerate(model(subvolumes.to(DEVICE))):\n                    output.append(value.tolist())\n\n        # Saving the files.\n        del eval_loader\n        np_output = np.array(output)\n        \n        print(\"The new mean is: \" + str(np.mean(np_output)))\n        print(\"The new std is: \" + str(np.std(np_output)))\n        \n        if np.std(np_output) > highest_avg:\n            highest_avg = np.std(np_output)\n            slices = image_selections\n        else:\n            print(\"A new best layer has not been found.\")\n\n        del np_output\n        gc.collect()\n        k += 2\n    print(\"The highest average is: \" +str(highest_avg))\n    print(\"The selected image layers are: \" + str(slices))\n    return slices","metadata":{"execution":{"iopub.status.busy":"2023-06-14T15:16:33.830640Z","iopub.execute_input":"2023-06-14T15:16:33.831856Z","iopub.status.idle":"2023-06-14T15:16:33.867422Z","shell.execute_reply.started":"2023-06-14T15:16:33.831808Z","shell.execute_reply":"2023-06-14T15:16:33.865958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This is detecting the locations of test fragments.\nbase_path = Path(\"/kaggle/input/vesuvius-challenge-ink-detection/test/\")\n#test_fragments = sorted([fragment_name for fragment_name in base_path.iterdir()])\ntest_fragments = [fragment_name for fragment_name in base_path.iterdir()]\n\n# This creates the output dictionary.\nsubmission_dict = {'ID': [], 'Predicted': []}","metadata":{"_uuid":"a264a491-7560-4abb-8364-2d94c6dda42c","_cell_guid":"6c6479bc-8062-4fff-b26f-ae5a8b2119c4","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-06-14T15:16:33.874456Z","iopub.execute_input":"2023-06-14T15:16:33.875547Z","iopub.status.idle":"2023-06-14T15:16:33.886399Z","shell.execute_reply.started":"2023-06-14T15:16:33.875503Z","shell.execute_reply":"2023-06-14T15:16:33.885168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for PREFIX in test_fragments:\n    fragment_name = PREFIX.name\n    PREFIX = str(PREFIX)+ '/'\n    print(PREFIX)\n\n    # This is processing the layers for the first pass.\n    image_selections = slice_finder(PREFIX)\n#     image_selections = [18,20,22,24,26,28,30,32,34,36]\n    img_layers = sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))\n\n    images = [np.array(Image.open(img_layers[i-1]), dtype=np.float32)/65535.0 for i in tqdm(image_selections)]\n    image_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n    del images\n\n    mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER_1:mask.shape[0]-BUFFER_1, BUFFER_1:mask.shape[1]-BUFFER_1] = True\n\n    grid = np.zeros(mask.shape)\n    i = WINDOW_1\n    while i < grid.shape[0]:\n        j = WINDOW_1\n        while j < grid.shape[1]:\n            grid[i,j] = 1\n            j += 2*WINDOW_1+1\n        i += 2*WINDOW_1+1\n\n    arr_mask = np.array(mask) * not_border * grid    \n    full_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    all_pixels = np.argwhere(full_rect)\n\n    del not_border, arr_mask, grid\n\n    eval_dataset = SubvolumeDataset(image_stack, all_pixels, BUFFER_1)\n\n    image_stack\n    gc.collect()\n\n    eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\n\n    # This is running the model on the first pass.\n    if torch.cuda.is_available():\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L1.pth').to(DEVICE)\n    else:\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L1.pth', map_location=torch.device('cpu')).to(DEVICE)\n    output_1=[]\n    model.eval()\n    with torch.no_grad():\n        for i, (subvolumes) in enumerate(eval_loader):\n            for j, value in enumerate(model(subvolumes.to(DEVICE))):\n                output_1.append(value.tolist())\n\n    # This is processing the layers for the second pass.\n    output_1 = np.array(output_1)\n\n    positives = 0\n    THRESHOLD = 1.00\n    while positives < 0.08:\n        THRESHOLD -= 0.01\n        threshold_test = np.where(output_1 > THRESHOLD, 1, 0)\n        positives = np.mean(threshold_test)\n\n    output_1 = torch.tensor(output_1).gt(THRESHOLD).cpu()\n    \n    layer_2_input = full_rect\n    for i in range(len(output_1)):\n        pixel_y = all_pixels[i][0]\n        pixel_x = all_pixels[i][1]\n        layer_2_input[pixel_y-WINDOW_1:pixel_y+WINDOW_1+1,pixel_x-WINDOW_1:pixel_x+WINDOW_1+1] = output_1[i].reshape(WINDOW_1*2+1,WINDOW_1*2+1)\n\n    # This is checking to make sure it looks normal.\n    fig, (ax1) = plt.subplots(1, 1)\n    ax1.imshow(layer_2_input, cmap='gray')\n    plt.show()\n\n    # This adds it to the stack.\n    image_stack = torch.concatenate((torch.tensor(layer_2_input, dtype=torch.float32).reshape(1,mask.shape[0],mask.shape[1]).to(DEVICE), image_stack))\n\n    # This defines the pixels we need to review based on the new buffer\n    mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER_2:mask.shape[0]-BUFFER_2, BUFFER_2:mask.shape[1]-BUFFER_2] = True\n\n    grid = np.zeros(mask.shape)\n    i = WINDOW_2\n    while i < grid.shape[0]:\n        j = WINDOW_2\n        while j < grid.shape[1]:\n            grid[i,j] = 1\n            j += 2*WINDOW_2+1\n        i += 2*WINDOW_2+1\n\n    arr_mask = np.array(mask) * not_border * grid    \n    full_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    all_pixels = np.argwhere(full_rect)\n    \n    del not_border, arr_mask, grid\n    \n    eval_dataset = SubvolumeDataset(image_stack, all_pixels, BUFFER_2)\n    \n    image_stack\n    gc.collect()\n    \n    eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)\n\n    # This is running the model on the second pass.\n    if torch.cuda.is_available():\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L2.pth').to(DEVICE)\n    else:\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L2.pth', map_location=torch.device('cpu')).to(DEVICE)\n    output_2=[]\n    model.eval()\n    with torch.no_grad():\n        for i, (subvolumes) in enumerate(eval_loader):\n            for j, value in enumerate(model(subvolumes.to(DEVICE))):\n                output_2.append(value.tolist())\n\n    # This is processing the new layers for the third pass.\n    output_2 = np.array(output_2)\n    \n    positives = 0\n    THRESHOLD = 1.00\n    while positives < 0.08:\n        THRESHOLD -= 0.01\n        threshold_test = np.where(output_2 > THRESHOLD, 1, 0)\n        positives = np.mean(threshold_test)\n\n    output_2 = torch.tensor(output_2).gt(THRESHOLD).cpu()\n    \n    layer_3_input = full_rect\n    for i in range(len(output_2)):\n        pixel_y = all_pixels[i][0]\n        pixel_x = all_pixels[i][1]\n        layer_3_input[pixel_y-WINDOW_2:pixel_y+WINDOW_2+1,pixel_x-WINDOW_2:pixel_x+WINDOW_2+1] = output_2[i].reshape(WINDOW_2*2+1,WINDOW_2*2+1)\n\n    # This is checking to make sure it looks normal.\n    fig, (ax1) = plt.subplots(1, 1)\n    ax1.imshow(layer_3_input, cmap='gray')\n    plt.show()\n                \n    # This is re-processing the layers for the final pass.\n        \n    image_selections = image_selections[1:len(image_selections)-1]\n    img_layers = sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))\n\n    images = [np.array(Image.open(img_layers[i-1]), dtype=np.float32)/65535.0 for i in tqdm(image_selections)]\n    image_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n    del images\n\n    mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER_3:mask.shape[0]-BUFFER_3, BUFFER_3:mask.shape[1]-BUFFER_3] = True\n\n    grid = np.zeros(mask.shape)\n    i = WINDOW_3\n    while i < grid.shape[0]:\n        j = WINDOW_3\n        while j < grid.shape[1]:\n            grid[i,j] = 1\n            j += 2*WINDOW_3+1\n        i += 2*WINDOW_3+1\n\n    arr_mask = np.array(mask) * not_border * grid    \n    full_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    all_pixels = np.argwhere(full_rect)\n\n    del not_border, arr_mask, grid  \n    \n    # This adds layer 2 to the stack.\n    image_stack = torch.concatenate((torch.tensor(layer_2_input, dtype=torch.float32).reshape(1,mask.shape[0],mask.shape[1]).to(DEVICE), image_stack))\n    # This adds layer 3 to the stack.\n    image_stack = torch.concatenate((torch.tensor(layer_3_input, dtype=torch.float32).reshape(1,mask.shape[0],mask.shape[1]).to(DEVICE), image_stack))\n    \n    eval_dataset = SubvolumeDataset(image_stack, all_pixels, BUFFER_3)\n    eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False) \n    \n    # This is running the model on the third pass.\n    if torch.cuda.is_available():\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L3.pth').to(DEVICE)\n    else:\n        model = torch.load('/kaggle/input/vevuvius-v5/v5_model_L3.pth', map_location=torch.device('cpu')).to(DEVICE)\n    output_3=[]\n    model.eval()\n    with torch.no_grad():\n        for i, (subvolumes) in enumerate(eval_loader):\n            for j, value in enumerate(model(subvolumes.to(DEVICE))):\n                output_3.append(value.tolist())\n    \n    # Here we're converting the output into a submittable file.\n    mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n    not_border = np.zeros(mask.shape, dtype=bool)\n    not_border[BUFFER_3:mask.shape[0]-BUFFER_3, BUFFER_3:mask.shape[1]-BUFFER_3] = True\n\n    grid = np.zeros(mask.shape)\n    i = WINDOW_3\n    while i < grid.shape[0]:\n        j = WINDOW_3\n        while j < grid.shape[1]:\n            grid[i,j] = 1\n            j += 2*WINDOW_3+1\n        i += 2*WINDOW_3+1\n\n    arr_mask = np.array(mask) * not_border * grid    \n    full_rect = np.ones(mask.shape, dtype=bool) * arr_mask\n    \n    \n    output_3 = np.array(output_3)\n    layer_4_input = full_rect\n\n    for i in range(len(output_3)):\n        pixel_y = all_pixels[i][0]\n        pixel_x = all_pixels[i][1]\n        layer_4_input[pixel_y-WINDOW_3:pixel_y+WINDOW_3+1,pixel_x-WINDOW_3:pixel_x+WINDOW_3+1] = output_3[i].reshape(WINDOW_3*2+1,WINDOW_3*2+1)\n\n    positives = 0\n    THRESHOLD = 1.00\n    while positives < 0.18:\n        THRESHOLD -= 0.01\n        threshold_test = np.where(layer_4_input > THRESHOLD, 1, 0)\n        positives = np.mean(threshold_test)\n#     THRESHOLD = .3 \n    print(\"The threshold is: \" + str(THRESHOLD))\n    \n    layer_4_input = torch.tensor(layer_4_input).gt(THRESHOLD).cpu()\n\n    # This is checking to make sure it looks normal.\n    fig, (ax1) = plt.subplots(1, 1)\n    ax1.imshow(layer_4_input, cmap='gray')\n    plt.show()\n\n    # This is converting the output into an RLE:\n    rle_output = rle(layer_4_input)\n    submission_dict['ID'].append(fragment_name)\n    submission_dict['Predicted'].append(rle_output)","metadata":{"_uuid":"eadd89c9-2b8a-4866-b17b-3258b3c61a21","_cell_guid":"ef66db80-a98a-4476-836d-155cac2e2caa","jupyter":{"outputs_hidden":false},"collapsed":false,"execution":{"iopub.status.busy":"2023-06-14T15:16:33.891352Z","iopub.execute_input":"2023-06-14T15:16:33.894197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(data=submission_dict)\nsubmission.to_csv('submission.csv', index=False)","metadata":{"_uuid":"5ffd5834-5955-464d-8564-e9f0ba4b5491","_cell_guid":"75ec4a94-671b-46c9-a398-1700471101d4","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(submission)","metadata":{"_uuid":"b14ac0a8-6db4-4aad-a28f-1ed862252dd2","_cell_guid":"b482321c-e965-45ff-9c5c-cd092fe6f44d","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n# !!! Everything below this line is testing, should be commented out. !!!\n# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!","metadata":{"_uuid":"9f926eff-b1d3-48cf-b198-8096294de7d4","_cell_guid":"b3525ffa-af8d-4d3c-9e4e-30b3cb48aab6","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]}]}