{"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":"# Surface Map\n\nThis notebook uses slightly different and I believe more accurate surface detection algorithm than by computing the \"minimum of the mean\" as proposed in original kernel: https://www.kaggle.com/code/raki21/surface-map-visualization\n\nInstead of calculating minimum of means across slices, sum of absolute differences between local mean and global mean is minimized. In addition, instead of calculating mean on tiles, which has edge effects even when calculating tiles with overlap, wide Gausian filter is used with sigma=200 (which corresponds to radius=800 and kernel=1601 https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html)\n\nThis was the main reason I precomputed Gausian filters for each layer of the train data. More on how it was done you can find below:\nhttps://www.kaggle.com/code/elvenmonk/build-light-maps-1-3\nhttps://www.kaggle.com/code/elvenmonk/build-light-maps-2","metadata":{}},{"cell_type":"markdown","source":"## Building Surface Maps","metadata":{}},{"cell_type":"code","source":"import os\nimport gc\nimport glob\nimport numpy as np\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nimport tifffile\n\nIN_PREFIX_1 = '/kaggle/input/build-light-maps-1-3/train/'\nIN_PREFIX_2 = '/kaggle/input/build-light-maps-2/train/'\nOUT_PREFIX = '/kaggle/working/train/'","metadata":{"execution":{"iopub.status.busy":"2023-05-30T06:35:43.324516Z","iopub.execute_input":"2023-05-30T06:35:43.326003Z","iopub.status.idle":"2023-05-30T06:35:43.597756Z","shell.execute_reply.started":"2023-05-30T06:35:43.325943Z","shell.execute_reply":"2023-05-30T06:35:43.596755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_image(fragment_id, z):\n    \"\"\"This function reads an image from a file and creates a foreground mask.\"\"\"\n\n    # Read image layers and store in list\n    IN_PREFIX = IN_PREFIX_2 if fragment_id == 2 else IN_PREFIX_1\n    for i, filename in enumerate(sorted(glob.glob(f'{IN_PREFIX}{fragment_id}/light/*.tif'))):\n        if i != z:\n            continue\n        image = tifffile.imread(filename)\n\n    return image\n\ndef calculate_surface_map(fragment_id, s=(10,40)):\n    \"\"\"This function calculates the surface map of the image.\"\"\"\n\n    # global means for each layer\n    d = s[1] - s[0]\n    means = np.zeros(d)\n    for z in tqdm(range(s[0],s[1])):\n        image = read_image(fragment_id, z)\n        shape = image.shape\n        means[z-s[0]] = image[image > 0].mean()\n    # find z shift of each pixel that aligns best with nearby pixels\n    v_distance = np.zeros((56-d, *shape), dtype=np.float32)\n    for z in tqdm(range(55)):\n        image = read_image(fragment_id, z)\n        for dz in range(max(0, z + d - 55), min(z+1, d)):\n            v_distance[z-dz] += np.abs(image - means[dz])\n\n    # Compute final surface map\n    surface_map = np.argmin(v_distance, axis=0)\n    del v_distance\n    gc.collect()\n    gc.collect()\n    \n    print(surface_map.shape, surface_map[surface_map > 0].min(), surface_map.max())\n    return surface_map\n\ndef plot_surface_map(surface_map):\n    \"\"\"This function plots the surface map and a selected slice of the image.\"\"\"\n    \n    # Plot the surface map on the left subplot\n    plt.figure(num=1, clear=True, figsize=(12,12))\n    surface_map_plot = plt.imshow(surface_map, cmap='plasma')\n    \n    # Add a colorbar legend to the surface map plot\n    color_bar = plt.colorbar(surface_map_plot)\n    color_bar.set_label('Elevation Value')\n    plt.show()\n\ndef visualize_surface(fragments=[2,3,1]):\n    \"\"\"This function visualizes the surface maps for a list of fragments.\"\"\"\n    \n    # Create a figure with subplots for each fragment\n    num_fragments = len(fragments)\n    \n    # Process each fragment and add to subplots\n    for fragment_id in fragments:\n        surface_map = calculate_surface_map(fragment_id)\n        if not os.path.exists(f'{OUT_PREFIX}{fragment_id}'):\n            os.makedirs(f'{OUT_PREFIX}{fragment_id}')\n        Image.fromarray(surface_map.astype(np.uint8)).save(f'{OUT_PREFIX}{fragment_id}/surface_map.png')\n        plot_surface_map(surface_map)\n","metadata":{"tags":[],"execution":{"iopub.status.busy":"2023-05-30T06:35:43.603257Z","iopub.execute_input":"2023-05-30T06:35:43.603653Z","iopub.status.idle":"2023-05-30T06:35:43.623408Z","shell.execute_reply.started":"2023-05-30T06:35:43.603621Z","shell.execute_reply":"2023-05-30T06:35:43.621980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_surface(fragments=[2,3,1])","metadata":{"tags":[],"execution":{"iopub.status.busy":"2023-05-30T06:35:43.629697Z","iopub.execute_input":"2023-05-30T06:35:43.630180Z"},"trusted":true},"execution_count":null,"outputs":[]}]}