{"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 Visualization\n\nThis notebook presents a surface mapping visualization for individual fragments. The aim is to estimate the location of the 'surface' in each of the 3 given fragments along the x and y axes. Notably, only fragment 2 seems to have significant variance in surface location. The surface estimation is done by computing the minimum of the mean across each slice for tiles in a fragment.\n\nThe notebook is divided into two parts:\n\n1. **Generating and displaying the surface maps**: In this section, we utilize image processing techniques to calculate and visualize the surface map of the given fragments. \n\n2. **Justifying the method for generating the map (min of mean across slices)**: This section provides a rationale and justification for our method of generating surface maps. We will explore the typical pattern that the mean across slices follows for each area in a fragment.","metadata":{}},{"cell_type":"markdown","source":"**Point of discussion:** Visually the surface maps seems to align well with the tif slices, but I have little knowledge about 3d-xray scans and might have the wrong intuition.","metadata":{}},{"cell_type":"markdown","source":"## Part 1: Showing Surface Maps","metadata":{}},{"cell_type":"code","source":"# Import necessary libraries\nimport os\nimport cv2\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-05-18T12:51:14.457547Z","iopub.execute_input":"2023-05-18T12:51:14.457996Z","iopub.status.idle":"2023-05-18T12:51:14.617287Z","shell.execute_reply.started":"2023-05-18T12:51:14.457952Z","shell.execute_reply":"2023-05-18T12:51:14.616059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Path to the competition data\ncomp_dir_path = '/kaggle/input/'\ncomp_folder_name = 'vesuvius-challenge-ink-detection'\ncomp_dataset_path = os.path.join(comp_dir_path, comp_folder_name)\n\n# Define the tile size and stride for image processing, these determine the resolution of the surface map\ntile_size = 512  # Can't be too small, explanation in Part 2\nstride = tile_size // 2","metadata":{"tags":[],"execution":{"iopub.status.busy":"2023-05-18T12:51:14.623128Z","iopub.execute_input":"2023-05-18T12:51:14.623498Z","iopub.status.idle":"2023-05-18T12:51:14.629829Z","shell.execute_reply.started":"2023-05-18T12:51:14.623466Z","shell.execute_reply":"2023-05-18T12:51:14.628389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_image(fragment_id, start_slice=20, end_slice=50):\n    \"\"\"This function reads an image from a file and creates a foreground mask.\"\"\"\n    \n    # Placeholder for storing image slices\n    images = []\n    idxs = range(start_slice, end_slice)\n    \n    # Read image slices and store in list\n    for i in tqdm(idxs):\n        image = cv2.imread(os.path.join(comp_dataset_path, f\"train/{fragment_id}/surface_volume/{i:02}.tif\"), 0)\n        images.append(image)\n        \n    # Stack image slices to create a 3D image\n    images = np.stack(images, axis=2)\n    \n    # Read and process the foreground mask\n    foreground = cv2.imread(os.path.join(comp_dataset_path, f\"train/{fragment_id}/mask.png\"), 0) // 255\n    return images, foreground\n\ndef calculate_surface_map(image, foreground, start_slice=20):\n    \"\"\"This function calculates the surface map of the image.\"\"\"\n    \n    # Define the ranges for x and y coordinates\n    x_range = range(0, image.shape[1]-tile_size+1, stride)\n    y_range = range(0, image.shape[0]-tile_size+1, stride)\n\n    # Initialize surface and overlap maps\n    surface_map = np.zeros((image.shape[0], image.shape[1]))\n    overlap_map = np.zeros((image.shape[0], image.shape[1])) + 1e-5\n    \n    # Compute surface map\n    for y1 in tqdm(y_range):\n        for x1 in x_range:\n            y2 = y1 + tile_size\n            x2 = x1 + tile_size\n            \n            # Only consider tiles that mainly contain scroll\n            if np.sum(foreground[y1:y2, x1:x2]) > tile_size**2 // 2: \n                tile = image[y1:y2, x1:x2]\n\n                # Calculate mean for each slice and find the minimum mean, we have some overlap between tiles \n                slice_means = np.mean(tile, axis=(0, 1))\n                min_mean = np.argmin(slice_means)\n                surface_map[y1:y2, x1:x2] += min_mean + start_slice   \n                overlap_map[y1:y2, x1:x2] += np.ones((tile_size, tile_size))\n\n    # Compute final surface map\n    surface_map = surface_map / overlap_map\n    return surface_map\n\ndef plot_surface_map(surface_map, image, axes, fig, fragment_index, slice_imshow=35, start_slice=20, end_slice=50):\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    surface_map_plot = axes[fragment_index, 0].imshow(surface_map, vmin=start_slice, vmax=end_slice)\n    axes[fragment_index, 0].set_title('Surface Map')\n    \n    # Add a colorbar legend to the surface map plot\n    color_bar = fig.colorbar(surface_map_plot, ax=axes[fragment_index, 0])\n    color_bar.set_label('Elevation Value')\n\n    # Load and plot the image slice on the right subplot\n    axes[fragment_index, 1].imshow(image[:, :, slice_imshow-start_slice])\n    axes[fragment_index, 1].set_title(f'Slice {slice_imshow} of the fragment')\n\ndef visualize_surface(fragments=[1,2,3], start_slice=20, end_slice=50, slice_imshow=35):\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    fig, axes = plt.subplots(num_fragments, 2, figsize=(20, 10*num_fragments))\n    \n    # Process each fragment and add to subplots\n    for fragment_index, fragment_id in enumerate(fragments):\n        image, foreground = read_image(fragment_id, start_slice=start_slice, end_slice=end_slice)\n        surface_map = calculate_surface_map(image, foreground, start_slice=start_slice)\n        plot_surface_map(surface_map, image, axes, fig, fragment_index, slice_imshow=slice_imshow, start_slice=start_slice, end_slice=end_slice)\n    \n    # Adjust the spacing between subplots for better visibility\n    plt.subplots_adjust(wspace=0.05)\n\n    # Display the plots\n    plt.show()","metadata":{"tags":[],"execution":{"iopub.status.busy":"2023-05-18T12:51:14.632222Z","iopub.execute_input":"2023-05-18T12:51:14.632656Z","iopub.status.idle":"2023-05-18T12:51:14.655630Z","shell.execute_reply.started":"2023-05-18T12:51:14.632606Z","shell.execute_reply":"2023-05-18T12:51:14.654321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example usage\nvisualize_surface(fragments=[1,2,3], start_slice=25, end_slice=55, slice_imshow=35)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2023-05-18T12:51:14.658808Z","iopub.execute_input":"2023-05-18T12:51:14.659256Z","iopub.status.idle":"2023-05-18T12:56:11.312971Z","shell.execute_reply.started":"2023-05-18T12:51:14.659219Z","shell.execute_reply":"2023-05-18T12:56:11.311699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Part 2: Method for Surface Calculation\n\nIn this section, we describe the approach used for estimating the surface location within each fragment. The mean across\nslices follows a typical pattern for each area in a fragment:\n\n1. A period of relatively constant values\n2. A peak in intensity\n3. A sharp decline across the next slices\n4. A valley, followed by an increase in values again\n\nWe will first illustrate this pattern by examining the entire fragments 1, 2, and 3. Afterward, we will inspect some\nsubareas within the fragments.","metadata":{}},{"cell_type":"code","source":"def calculate_stats_for_tile(image, top_left_x, top_left_y, tile_size, start_slice=25):\n    \"\"\"This function calculates statistics for a specific tile and returns the mean values for each slice.\"\"\"\n\n    # Extract the tile from the image\n    if tile_size == 'full':\n        tile = image\n    else:\n        tile = image[top_left_x:top_left_x+tile_size, top_left_y:top_left_y+tile_size]  \n    \n    # Calculate mean for each slice\n    slice_means = np.mean(tile, axis=(0,1))\n\n    return slice_means\n\ndef plot_images_and_mean_slices(fragments=[1, 2, 3], start_slice=10, end_slice=60, tile_size='full', slice_imshow=35):\n    \"\"\"This function plots images and mean slice plots for the given fragments in a 2x3 figure.\"\"\"\n\n    # Create a 2x3 figure\n    fig, axs = plt.subplots(2, 3, figsize=(20, 10))\n\n    # For each fragment\n    for i, fragment in enumerate(fragments):\n        # Load the image\n        image, _ = read_image(fragment, start_slice=start_slice, end_slice=end_slice)\n\n        # Plot the slice number 35 of the image on the top row\n        axs[0, i].imshow(image[:, :, slice_imshow-start_slice])\n        axs[0, i].set_title(f'Fragment {fragment}: Slice 35')\n\n        # Calculate and plot the mean of slices on the bottom row\n        slice_means = calculate_stats_for_tile(image, top_left_x=0, top_left_y=0, tile_size=tile_size, start_slice=start_slice)\n        axs[1, i].plot(range(start_slice, start_slice+len(slice_means)), slice_means, marker='o')\n        axs[1, i].set_title(f'Fragment {fragment}: Mean of slices')\n        axs[1, i].set_xlabel('Slice')\n        axs[1, i].set_ylabel('Mean')\n\n    # Show the figure\n    plt.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-18T12:56:11.314524Z","iopub.execute_input":"2023-05-18T12:56:11.315610Z","iopub.status.idle":"2023-05-18T12:56:11.331663Z","shell.execute_reply.started":"2023-05-18T12:56:11.315561Z","shell.execute_reply":"2023-05-18T12:56:11.330472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Call the function to plot images and mean slice plots for fragments 1, 2 and 3\nplot_images_and_mean_slices(start_slice=10, end_slice=60, tile_size='full', slice_imshow=35)","metadata":{"execution":{"iopub.status.busy":"2023-05-18T12:56:11.333718Z","iopub.execute_input":"2023-05-18T12:56:11.334599Z","iopub.status.idle":"2023-05-18T13:02:04.533301Z","shell.execute_reply.started":"2023-05-18T12:56:11.334557Z","shell.execute_reply":"2023-05-18T13:02:04.531925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Subarea Showcase\n\nThe following figures demonstrate the variability of the mean curve based on the size of the subsection. You will notice that as the subsection decreases in size, the mean curve becomes increasingly chaotic. \n\nThis variability ultimately imposes a limit on the resolution achievable in the surface map.","metadata":{}},{"cell_type":"code","source":"def calculate_and_plot_stats(image, axs, top_left_x, top_left_y, tile_size, start_slice=25, slice_imshow=35):\n    \"\"\"This function calculates statistics for a specific tile and plots the mean values for each slice.\"\"\"\n\n    # Extract the tile from the image\n    if tile_size == 'full':\n        tile = image\n    else:\n        tile = image[top_left_x:top_left_x+tile_size, top_left_y:top_left_y+tile_size]  \n    \n    # Calculate mean for each slice\n    slice_means = np.mean(tile, axis=(0,1))\n\n    # Calculate slice with highest and lowest mean\n    max_mean_idx = np.argmax(slice_means)\n    min_mean_idx = np.argmin(slice_means)\n    print(f\"Slice with highest mean: {max_mean_idx + start_slice}\")\n    print(f\"Slice with lowest mean: {min_mean_idx + start_slice}\")\n\n    # Plot mean\n    axs[0].plot(range(start_slice, start_slice+len(slice_means)), slice_means, marker='o')\n    axs[0].scatter([max_mean_idx + start_slice, min_mean_idx + start_slice], [slice_means[max_mean_idx], slice_means[min_mean_idx]], color='red')  # mark highest and lowest\n    axs[0].set_title('Mean of each slice')\n    axs[0].set_xlabel('Slice')\n    axs[0].set_ylabel('Mean')\n\n    # Plot slice\n    axs[1].imshow(tile[:,:,slice_imshow-start_slice])\n    axs[1].set_title(f\"Slice {slice_imshow} of the fragment\")","metadata":{"execution":{"iopub.status.busy":"2023-05-18T13:02:04.535076Z","iopub.execute_input":"2023-05-18T13:02:04.536411Z","iopub.status.idle":"2023-05-18T13:02:04.551216Z","shell.execute_reply.started":"2023-05-18T13:02:04.536365Z","shell.execute_reply":"2023-05-18T13:02:04.549849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Analyze subsections of fragment 2 to observe how the curve of the mean for each slice changes with different resolutions\nstart_slice = 10\ntile_sizes = [4000, 2000, 1000, 500, 250, 20]\nimage_2, _ = read_image(2, start_slice=start_slice, end_slice=60)\n\n# Calculate the number of subplots needed (2 figures per tile_size)\nn_subplots = len(tile_sizes) * 2\n\n# Calculate the number of rows needed for the subplots (2 subplots per row)\nn_rows = len(tile_sizes)\n\n# Create a figure with the calculated number of subplots\nfig, axs = plt.subplots(n_rows, 2, figsize=(15, 5 * n_rows))\n\nfor idx, tile_size in enumerate(tile_sizes):\n    print(f'\\nAnalysing tile of size: {tile_size}\\n')\n    \n    # Calculate and plot stats, passing the current row of subplots\n    calculate_and_plot_stats(image_2, axs[idx], top_left_x=10000, top_left_y=3500, tile_size=tile_size, start_slice=start_slice, slice_imshow=35)\n\n# Adjust the spacing between subplots and show the plot\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-18T13:02:29.744531Z","iopub.execute_input":"2023-05-18T13:02:29.745002Z","iopub.status.idle":"2023-05-18T13:05:47.344256Z","shell.execute_reply.started":"2023-05-18T13:02:29.744967Z","shell.execute_reply":"2023-05-18T13:05:47.342574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}