{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "import open3d as o3d\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import os\n", "import sys\n", "\n", "# only needed for tutorial, monkey patches visualization\n", "sys.path.append('..')\n", "import open3d_tutorial as o3dtut\n", "# change to True if you want to interact with the visualization windows\n", "o3dtut.interactive = not \"CI\" in os.environ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Surface reconstruction\n", "\n", "In many scenarios we want to generate a dense 3D geometry, i.e., a triangle mesh. However, from a multi-view stereo method, or a depth sensor we only obtain an unstructured point cloud. To get a triangle mesh from this unstructured input we need to perform surface reconstruction. In the literature there exists a couple of methods and Open3D currently implements the following:\n", "\n", "- Alpha shapes [\\[Edelsbrunner1983\\]](../reference.html#Edelsbrunner1983)\n", "- Ball pivoting [\\[Bernardini1999\\]](../reference.html#Bernardini1999)\n", "- Poisson surface reconstruction [\\[Kazhdan2006\\]](../reference.html#Kazhdan2006)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Alpha shapes\n", "The alpha shape [\\[Edelsbrunner1983\\]](../reference.html#Edelsbrunner1983) is a generalization of a convex hull. As described [here](https://graphics.stanford.edu/courses/cs268-11-spring/handouts/AlphaShapes/as_fisher.pdf) one can intuitively\n", "think of an alpha shape as the following: Imagine a huge mass of ice cream containing the points $S$ as hard chocolate pieces. Using one of these sphere-formed ice cream spoons we carve out all parts of the ice cream block we can reach without bumping into chocolate pieces, thereby even carving out holes in the inside (e.g., parts not reachable by simply moving the\n", "spoon from the outside). We will eventually end up with a (not necessarily convex) object bounded by caps, arcs and points. If we now straighten all round faces to triangles and line segments, we have an intuitive description of what is called the alpha shape of $S$.\n", "\n", "Open3D implements the method `create_from_point_cloud_alpha_shape` that involves the tradeoff parameter `alpha`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bunny = o3d.data.BunnyMesh()\n", "mesh = o3d.io.read_triangle_mesh(bunny.path)\n", "mesh.compute_vertex_normals()\n", "\n", "pcd = mesh.sample_points_poisson_disk(750)\n", "o3d.visualization.draw_geometries([pcd])\n", "alpha = 0.03\n", "print(f\"alpha={alpha:.3f}\")\n", "mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha)\n", "mesh.compute_vertex_normals()\n", "o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The implementation is based on the convex hull of the point cloud. If we want to compute multiple alpha shapes from a given point cloud, then we can save some computation by only computing the convex hull once and pass it to `create_from_point_cloud_alpha_shape`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tetra_mesh, pt_map = o3d.geometry.TetraMesh.create_from_point_cloud(pcd)\n", "for alpha in np.logspace(np.log10(0.5), np.log10(0.01), num=4):\n", " print(f\"alpha={alpha:.3f}\")\n", " mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(\n", " pcd, alpha, tetra_mesh, pt_map)\n", " mesh.compute_vertex_normals()\n", " o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ball pivoting\n", "The ball pivoting algorithm (BPA) [\\[Bernardini1999\\]](../reference.html#Bernardini1999) is a surface reconstruction method which is related to alpha shapes. Intuitively, think of a 3D ball with a given radius that we drop on the point cloud. If it hits any 3 points (and it does not fall through those 3 points) it creates a triangles. Then, the algorithm starts pivoting from the edges of the existing triangles and every time it hits 3 points where the ball does not fall through we create another triangle.\n", "\n", "Open3D implements this method in `create_from_point_cloud_ball_pivoting`. The method accepts a list of `radii` as parameter that corresponds to the radii of the individual balls that are pivoted on the point cloud.\n", "\n", "\n", "
\n", " \n", "**Note:** \n", "\n", "This algorithm assumes that the `PointCloud` has `normals`.\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bunny = o3d.data.BunnyMesh()\n", "gt_mesh = o3d.io.read_triangle_mesh(bunny.path)\n", "gt_mesh.compute_vertex_normals()\n", "\n", "pcd = gt_mesh.sample_points_poisson_disk(3000)\n", "o3d.visualization.draw_geometries([pcd])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "radii = [0.005, 0.01, 0.02, 0.04]\n", "rec_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(\n", " pcd, o3d.utility.DoubleVector(radii))\n", "o3d.visualization.draw_geometries([pcd, rec_mesh])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Poisson surface reconstruction\n", "The Poisson surface reconstruction method [\\[Kazhdan2006\\]](../reference.html#Kazhdan2006) solves a regularized optimization problem to obtain a smooth surface. For this reason, Poisson surface reconstruction can be preferable to the methods mentioned above, as they produce non-smooth results since the points of the `PointCloud` are also the `vertices` of the resulting triangle mesh without any modifications.\n", "\n", "Open3D implements the method `create_from_point_cloud_poisson` which is basically a wrapper of the code of [Kazhdan](https://github.com/mkazhdan/PoissonRecon). An important parameter of the function is `depth` that defines the depth of the octree used for the surface reconstruction and hence implies the resolution of the resulting triangle mesh. A higher `depth` value means a mesh with more details.\n", "\n", "
\n", " \n", "**Note:** \n", "\n", "This algorithm assumes that the `PointCloud` has `normals`.\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eagle = o3d.data.EaglePointCloud()\n", "pcd = o3d.io.read_point_cloud(eagle.path)\n", "\n", "print(pcd)\n", "o3d.visualization.draw_geometries([pcd],\n", " zoom=0.664,\n", " front=[-0.4761, -0.4698, -0.7434],\n", " lookat=[1.8900, 3.2596, 0.9284],\n", " up=[0.2304, -0.8825, 0.4101])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('run Poisson surface reconstruction')\n", "with o3d.utility.VerbosityContextManager(\n", " o3d.utility.VerbosityLevel.Debug) as cm:\n", " mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(\n", " pcd, depth=9)\n", "print(mesh)\n", "o3d.visualization.draw_geometries([mesh],\n", " zoom=0.664,\n", " front=[-0.4761, -0.4698, -0.7434],\n", " lookat=[1.8900, 3.2596, 0.9284],\n", " up=[0.2304, -0.8825, 0.4101])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Poisson surface reconstruction will also create triangles in areas of low point density, and even extrapolates into some areas (see bottom of the eagle output above). The `create_from_point_cloud_poisson` function has a second `densities` return value that indicates for each vertex the density. A low density value means that the vertex is only supported by a low number of points from the input point cloud.\n", "\n", "In the code below we visualize the density in 3D using pseudo color. Violet indicates low density and yellow indicates a high density." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('visualize densities')\n", "densities = np.asarray(densities)\n", "density_colors = plt.get_cmap('plasma')(\n", " (densities - densities.min()) / (densities.max() - densities.min()))\n", "density_colors = density_colors[:, :3]\n", "density_mesh = o3d.geometry.TriangleMesh()\n", "density_mesh.vertices = mesh.vertices\n", "density_mesh.triangles = mesh.triangles\n", "density_mesh.triangle_normals = mesh.triangle_normals\n", "density_mesh.vertex_colors = o3d.utility.Vector3dVector(density_colors)\n", "o3d.visualization.draw_geometries([density_mesh],\n", " zoom=0.664,\n", " front=[-0.4761, -0.4698, -0.7434],\n", " lookat=[1.8900, 3.2596, 0.9284],\n", " up=[0.2304, -0.8825, 0.4101])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can further use the density values to remove vertices and triangles that have a low support. In the code below we remove all vertices (and connected triangles) that have a lower density value than the $0.01$ quantile of all density values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('remove low density vertices')\n", "vertices_to_remove = densities < np.quantile(densities, 0.01)\n", "mesh.remove_vertices_by_mask(vertices_to_remove)\n", "print(mesh)\n", "o3d.visualization.draw_geometries([mesh],\n", " zoom=0.664,\n", " front=[-0.4761, -0.4698, -0.7434],\n", " lookat=[1.8900, 3.2596, 0.9284],\n", " up=[0.2304, -0.8825, 0.4101])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Normal estimation\n", "In the examples above we assumed that the point cloud has normals that point outwards. However, not all point clouds already come with associated normals. Open3D can be used to estimate point cloud normals with `estimate_normals`, which locally fits a plane per 3D point to derive the normal. However, the estimated normals might not be consistently oriented. `orient_normals_consistent_tangent_plane` propagates the normal orientation using a minimum spanning tree." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bunny = o3d.data.BunnyMesh()\n", "gt_mesh = o3d.io.read_triangle_mesh(bunny.path)\n", "\n", "pcd = gt_mesh.sample_points_poisson_disk(5000)\n", "pcd.normals = o3d.utility.Vector3dVector(np.zeros(\n", " (1, 3))) # invalidate existing normals\n", "\n", "pcd.estimate_normals()\n", "o3d.visualization.draw_geometries([pcd], point_show_normal=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pcd.orient_normals_consistent_tangent_plane(100)\n", "o3d.visualization.draw_geometries([pcd], point_show_normal=True)" ] } ], "metadata": { "celltoolbar": "Edit Metadata", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }