{"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":"# **Point Cloud conversion**\n\nThis notebook shows how to convert the voxels into a Point Cloud. A point cloud is an unordered set able to represent 3D structure, while a voxel represents a value on a 3D grid and can be akin to a 3D pixel. \n\nModeling the data as a point cloud as opposed to voxels could allow for structural features to be exploited in a different manner. [](http://)","metadata":{}},{"cell_type":"code","source":"import os\nfrom glob import glob\nfrom tifffile import tifffile\nimport numpy as np\nimport PIL.Image as Image\nimport torch.utils.data as data\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed","metadata":{"execution":{"iopub.status.busy":"2023-04-15T12:57:58.994504Z","iopub.execute_input":"2023-04-15T12:57:58.995711Z","iopub.status.idle":"2023-04-15T12:58:01.212776Z","shell.execute_reply.started":"2023-04-15T12:57:58.995648Z","shell.execute_reply":"2023-04-15T12:58:01.211318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Install libraries","metadata":{}},{"cell_type":"code","source":"!pip install open3d","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import open3d as o3","metadata":{"execution":{"iopub.status.busy":"2023-04-15T13:16:53.275177Z","iopub.execute_input":"2023-04-15T13:16:53.275564Z","iopub.status.idle":"2023-04-15T13:16:53.280477Z","shell.execute_reply.started":"2023-04-15T13:16:53.275530Z","shell.execute_reply":"2023-04-15T13:16:53.279757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/1/'\nBUFFER = 30  # Buffer size in x and y direction\nZ_START = 27 # First slice in the z direction to use\nZ_DIM = 10   # Number of slices in the z direction","metadata":{"execution":{"iopub.status.busy":"2023-04-15T12:58:11.457599Z","iopub.execute_input":"2023-04-15T12:58:11.458135Z","iopub.status.idle":"2023-04-15T12:58:11.463370Z","shell.execute_reply.started":"2023-04-15T12:58:11.458103Z","shell.execute_reply":"2023-04-15T12:58:11.462249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"mask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\nplt.imshow(mask, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2023-04-15T14:50:05.868467Z","iopub.execute_input":"2023-04-15T14:50:05.868894Z","iopub.status.idle":"2023-04-15T14:50:07.921413Z","shell.execute_reply.started":"2023-04-15T14:50:05.868855Z","shell.execute_reply":"2023-04-15T14:50:07.920457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the 3d x-ray scan, one slice at a time\nimages = [(tifffile.imread(filename)/65535.0).astype(np.float32) for filename \n          in tqdm(sorted(glob(PREFIX+\"surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM])]","metadata":{"execution":{"iopub.status.busy":"2023-04-15T13:13:48.804232Z","iopub.execute_input":"2023-04-15T13:13:48.804642Z","iopub.status.idle":"2023-04-15T13:13:50.455019Z","shell.execute_reply.started":"2023-04-15T13:13:48.804606Z","shell.execute_reply":"2023-04-15T13:13:50.453910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Point Cloud\n\nTo create the point cloud, we will unformly sample from each masked surface volume. The point cloud mainly provides structural information in 3D space. We also have a single native feature of the point cloud, and that is the intensity of each point. For this notebook, we will just use the intenisty as the display color.\n\n\n","metadata":{}},{"cell_type":"code","source":"ROWS = images[0].shape[0]\nCOLS = images[0].shape[1]\nN_SAMPLES = 10000","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:35:47.978397Z","iopub.execute_input":"2023-04-15T15:35:47.978805Z","iopub.status.idle":"2023-04-15T15:35:47.984457Z","shell.execute_reply.started":"2023-04-15T15:35:47.978768Z","shell.execute_reply":"2023-04-15T15:35:47.983252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample from each surface volume\n\nInspired from: https://stackoverflow.com/questions/60384541/drawing-a-random-sample-from-a-numpy-array-with-index","metadata":{}},{"cell_type":"code","source":"# sample from valid regions of surface volume\nc = mask.ravel().cumsum()\nsamples = np.random.uniform(low=0, \n                            high=c[-1], \n                            size=(N_SAMPLES, Z_DIM)).astype(int)\n# get valid indexes\nx, y = np.unravel_index(c.searchsorted(samples), mask.shape)\nx, y = x[np.newaxis, ...], y[np.newaxis, ...]\n\n# get z dimensions from surface volume locations\nz = np.arange(Z_START, Z_START+Z_DIM)\nz = np.tile(z, N_SAMPLES).reshape(N_SAMPLES, -1)[np.newaxis, ...]","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:35:48.566464Z","iopub.execute_input":"2023-04-15T15:35:48.566851Z","iopub.status.idle":"2023-04-15T15:35:48.934315Z","shell.execute_reply.started":"2023-04-15T15:35:48.566818Z","shell.execute_reply":"2023-04-15T15:35:48.933095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get point cloud\nxyz = np.vstack((x, y, z))\nxyz.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:35:48.936162Z","iopub.execute_input":"2023-04-15T15:35:48.936859Z","iopub.status.idle":"2023-04-15T15:35:48.946016Z","shell.execute_reply.started":"2023-04-15T15:35:48.936813Z","shell.execute_reply":"2023-04-15T15:35:48.944420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample from surface volumes to get intensities","metadata":{}},{"cell_type":"code","source":"intensities = np.zeros((N_SAMPLES, Z_DIM))\nfor i, img in enumerate(images):\n    intensities[:, i] = img[xyz[0, :, i], xyz[1, :, i]]","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:35:49.576753Z","iopub.execute_input":"2023-04-15T15:35:49.577151Z","iopub.status.idle":"2023-04-15T15:35:49.588009Z","shell.execute_reply.started":"2023-04-15T15:35:49.577116Z","shell.execute_reply":"2023-04-15T15:35:49.586546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Need to reshape. Optional: check that everything lines up across surface volumes","metadata":{}},{"cell_type":"code","source":"print(xyz[:, 20, 1], intensities[20, 1])\nprint(xyz.T.reshape((-1, 3))[20 + N_SAMPLES, :], intensities.T.reshape((-1))[20 + N_SAMPLES])","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:22:35.016462Z","iopub.execute_input":"2023-04-15T15:22:35.016869Z","iopub.status.idle":"2023-04-15T15:22:35.024525Z","shell.execute_reply.started":"2023-04-15T15:22:35.016835Z","shell.execute_reply":"2023-04-15T15:22:35.023188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Reshape and Normalize","metadata":{}},{"cell_type":"code","source":"xyz = xyz.T.reshape((-1, 3))\nxyz = xyz/xyz.max(axis=0)\n\nintensities = intensities.T.reshape((-1)).repeat((3)).reshape((-1, 3))","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:22:36.215845Z","iopub.execute_input":"2023-04-15T15:22:36.216180Z","iopub.status.idle":"2023-04-15T15:22:36.228252Z","shell.execute_reply.started":"2023-04-15T15:22:36.216151Z","shell.execute_reply":"2023-04-15T15:22:36.226780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convert to open3d point cloud object\n\nWe will use a colormap to enhance the viz","metadata":{}},{"cell_type":"code","source":"colors = plt.get_cmap('cool')\ncolors","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:29:20.251995Z","iopub.execute_input":"2023-04-15T15:29:20.252364Z","iopub.status.idle":"2023-04-15T15:29:20.266963Z","shell.execute_reply.started":"2023-04-15T15:29:20.252331Z","shell.execute_reply":"2023-04-15T15:29:20.265729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pcd = o3.geometry.PointCloud()\npcd.points = o3.utility.Vector3dVector(xyz) \npcd.colors = o3.utility.Vector3dVector(colors(intensities)[:, 0, :3])\n\npcd","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:29:20.766656Z","iopub.execute_input":"2023-04-15T15:29:20.767256Z","iopub.status.idle":"2023-04-15T15:29:20.831264Z","shell.execute_reply.started":"2023-04-15T15:29:20.767214Z","shell.execute_reply":"2023-04-15T15:29:20.830291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Display Point Cloud","metadata":{}},{"cell_type":"code","source":"o3.visualization.draw_plotly([pcd])","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:29:21.521731Z","iopub.execute_input":"2023-04-15T15:29:21.522087Z","iopub.status.idle":"2023-04-15T15:29:21.629523Z","shell.execute_reply.started":"2023-04-15T15:29:21.522054Z","shell.execute_reply":"2023-04-15T15:29:21.628273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(1, len(images), figsize=(15, 3))\nfor image, ax in zip(images, axes):\n    ax.imshow(colors(np.array(Image.fromarray(image).resize((image.shape[1]//20, image.shape[0]//20)), dtype=np.float32)))\n    ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-15T15:34:20.904052Z","iopub.execute_input":"2023-04-15T15:34:20.905031Z","iopub.status.idle":"2023-04-15T15:34:23.955536Z","shell.execute_reply.started":"2023-04-15T15:34:20.904979Z","shell.execute_reply":"2023-04-15T15:34:23.954435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}