{"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":"markdown","source":"# Preamble\n\nWhen looking at raw CT scan images it feels like difference between nearly straight lines of papyrus texture, and curly spots in many places where ink is labelled can be a key to solving this puzzle.\n\nThis gave me an idea to treat signal intensity as if it was elevation in DEM GIS model and calculate slope and more importantly slope direction (aspect) in each pixel to see if anything interesting can be seen there.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport glob\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom io import StringIO\nfrom sklearn.metrics import fbeta_score\nfrom skimage.util import view_as_windows\nfrom scipy.ndimage import distance_transform_edt, gaussian_filter\nfrom numba import jit","metadata":{"execution":{"iopub.status.busy":"2023-06-05T07:36:00.053482Z","iopub.execute_input":"2023-06-05T07:36:00.054962Z","iopub.status.idle":"2023-06-05T07:36:00.063627Z","shell.execute_reply.started":"2023-06-05T07:36:00.054905Z","shell.execute_reply":"2023-06-05T07:36:00.061710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Constants\nPREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/3/'\n\n# Load mask image\nmask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n# Load label image\nlabel = (np.array(Image.open(PREFIX+\"inklabels.png\")) > 0).astype(np.float32)\n\n# Load infrared image\nir = np.array(Image.open(PREFIX+\"ir.png\")).astype(np.float32)","metadata":{"execution":{"iopub.status.busy":"2023-06-05T07:25:22.322203Z","iopub.execute_input":"2023-06-05T07:25:22.322673Z","iopub.status.idle":"2023-06-05T07:25:23.217880Z","shell.execute_reply.started":"2023-06-05T07:25:22.322639Z","shell.execute_reply":"2023-06-05T07:25:23.216524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install xarray-spatial","metadata":{"execution":{"iopub.status.busy":"2023-06-05T07:08:16.650647Z","iopub.execute_input":"2023-06-05T07:08:16.651042Z","iopub.status.idle":"2023-06-05T07:08:34.613500Z","shell.execute_reply.started":"2023-06-05T07:08:16.651013Z","shell.execute_reply":"2023-06-05T07:08:34.612442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import xarray as xr\nfrom xrspatial.aspect import aspect\nfrom xrspatial.slope import slope\n\ndef plot_horizontal(images):\n    fig, ax = plt.subplots(1, len(images), figsize=(14,5))\n    for i, image in enumerate(images):\n        ax[i].imshow(image[4000:6000,2000:4000], cmap='binary')\n    plt.show()\n    \ndef plot_slope_aspect_maps(image):\n    xarr = xr.DataArray(np.pad(gaussian_filter(image, 10), 1))\n    image_asp = aspect(xarr)[1:-1,1:-1]/360.\n    image_slp = slope(xarr)[1:-1,1:-1]/90.\n    print(image_asp.data.min(),image_asp.data.max(),image_slp.data.min(),image_slp.data.max())\n    plot_horizontal([image, image_asp.data, image_slp.data, label])\n    \nplot_slope_aspect_maps(ir)","metadata":{"execution":{"iopub.status.busy":"2023-06-05T07:37:06.357373Z","iopub.execute_input":"2023-06-05T07:37:06.357840Z","iopub.status.idle":"2023-06-05T07:37:16.497895Z","shell.execute_reply.started":"2023-06-05T07:37:06.357805Z","shell.execute_reply":"2023-06-05T07:37:16.496542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Interesting, right? Slope and aspect maps still preserve a lot of ink information, even though ir image does not represent geometrical features.","metadata":{}},{"cell_type":"code","source":"for z, filename in tqdm(enumerate(sorted(glob.glob(PREFIX+\"surface_volume/*.tif\")))):\n    img = np.array(Image.open(filename), dtype=np.float32) / 257.\n    plot_slope_aspect_maps(img)","metadata":{"execution":{"iopub.status.busy":"2023-06-05T07:37:54.524007Z","iopub.execute_input":"2023-06-05T07:37:54.524446Z","iopub.status.idle":"2023-06-05T07:38:23.730862Z","shell.execute_reply.started":"2023-06-05T07:37:54.524398Z","shell.execute_reply":"2023-06-05T07:38:23.728646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Can you see anything?","metadata":{}},{"cell_type":"code","source":"# Function to generate run-length encoding (RLE) for the binary mask\ndef rle(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\n\n# Generate RLE for the binary output\nrle_output = rle(mask)\n\n# Save the RLE to a CSV file for submission\nwith open('submission.csv', 'w') as f:\n    f.write(\"Id,Predicted\\n\")\n    f.write(\"a,\" + rle_output + \"\\n\")\n    f.write(\"b,\" + rle_output + \"\\n\")\n\nprint(\"Submission file 'submission.csv' has been generated.\")","metadata":{"execution":{"iopub.status.busy":"2023-05-17T05:26:34.407695Z","iopub.status.idle":"2023-05-17T05:26:34.408326Z","shell.execute_reply.started":"2023-05-17T05:26:34.408033Z","shell.execute_reply":"2023-05-17T05:26:34.408061Z"},"trusted":true},"execution_count":null,"outputs":[]}]}