{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import open3d as o3d\n",
    "import numpy as np\n",
    "import copy\n",
    "import os\n",
    "import sys\n",
    "\n",
    "# monkey patches visualization and provides helpers to load geometries\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": [
    "# Mesh\n",
    "Open3D has a data structure for 3D triangle meshes called `TriangleMesh`.\n",
    "The code below shows how to read a triangle mesh from a `ply` and prints its vertices and triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Testing mesh in open3d ...\")\n",
    "mesh = o3dtut.get_knot_mesh()\n",
    "print(mesh)\n",
    "print('Vertices:')\n",
    "print(np.asarray(mesh.vertices))\n",
    "print('Triangles:')\n",
    "print(np.asarray(mesh.triangles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `TriangleMesh` class has a few data fields such as `vertices` and `triangles`. Open3D provides direct memory access to these fields via numpy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize a 3D mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Try to render a mesh with normals (exist: \" +\n",
    "        str(mesh.has_vertex_normals()) + \") and colors (exist: \" +\n",
    "        str(mesh.has_vertex_colors()) + \")\")\n",
    "o3d.visualization.draw_geometries([mesh])\n",
    "print(\"A mesh with no normals and no colors does not look good.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can rotate and move the mesh but it is painted with uniform gray color and does not look “3d”. The reason is that the current mesh does not have normals for vertices or faces. So uniform color shading is used instead of a more sophisticated Phong shading."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Surface normal estimation\n",
    "Let’s draw the mesh with surface normals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Computing normal and rendering it.\")\n",
    "mesh.compute_vertex_normals()\n",
    "print(np.asarray(mesh.triangle_normals))\n",
    "o3d.visualization.draw_geometries([mesh])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It uses `compute_vertex_normals` and `paint_uniform_color` which are member functions of `mesh`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crop mesh\n",
    "We remove half of the surface by directly operating on the t`riangle` and `triangle_normals` data fields of the mesh. This is done via numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"We make a partial mesh of only the first half triangles.\")\n",
    "mesh1 = copy.deepcopy(mesh)\n",
    "mesh1.triangles = o3d.utility.Vector3iVector(\n",
    "        np.asarray(mesh1.triangles)[:len(mesh1.triangles) // 2, :])\n",
    "mesh1.triangle_normals = o3d.utility.Vector3dVector(\n",
    "    np.asarray(mesh1.triangle_normals)[:len(mesh1.triangle_normals) // 2, :])\n",
    "print(mesh1.triangles)\n",
    "o3d.visualization.draw_geometries([mesh1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paint mesh\n",
    "Painting mesh is the same as how it worked for point cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Painting the mesh\")\n",
    "mesh1.paint_uniform_color([1, 0.706, 0])\n",
    "o3d.visualization.draw_geometries([mesh1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh properties\n",
    "A triangle mesh has several properties that can be tested with Open3D. One important property is the manifold property, where we can test the triangle mesh if it is edge manifold `is_edge_manifold` and if it is `is_vertex_manifold`. A triangle mesh is edge manifold, if each edge is bounding either one or two triangles. The function `is_edge_manifold` has the `bool` parameter `allow_boundary_edges` that defines if boundary edges should be allowed. Further, a triangle mesh is vertex manifold if the star of the vertex is edge‐manifold and edge‐connected, e.g., two or more faces connected only by a vertex and not by an edge.\n",
    "\n",
    "Another property is the test of self-intersection. The function `is_self_intersecting` returns `True` if there exists a triangle in the mesh that is intersecting another mesh. A watertight mesh can be defined as a mesh that is edge manifold, vertex manifold and not self intersecting. The function `is_watertight` implements this check in Open3D.\n",
    "\n",
    "We also can test the triangle mesh, if it is orientable, i.e. the triangles can be oriented in such a way that all normals point towards the outside. The corresponding function in Open3D is called `is_orientable`.\n",
    "\n",
    "The code below test a number of triangle meshes against those properties and visualizes the results. Non-manifold edges are shown in red, boundary edges in green, non-manifold vertices are visualized as green points, and self-intersecting triangles are shown in pink. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_properties(name, mesh):\n",
    "    mesh.compute_vertex_normals()\n",
    "    \n",
    "    edge_manifold = mesh.is_edge_manifold(allow_boundary_edges=True)\n",
    "    edge_manifold_boundary = mesh.is_edge_manifold(allow_boundary_edges=False)\n",
    "    vertex_manifold = mesh.is_vertex_manifold()\n",
    "    self_intersecting = mesh.is_self_intersecting()\n",
    "    watertight = mesh.is_watertight()\n",
    "    orientable = mesh.is_orientable()\n",
    "    \n",
    "    print(name)\n",
    "    print(f\"  edge_manifold:          {edge_manifold}\")\n",
    "    print(f\"  edge_manifold_boundary: {edge_manifold_boundary}\")\n",
    "    print(f\"  vertex_manifold:        {vertex_manifold}\")\n",
    "    print(f\"  self_intersecting:      {self_intersecting}\")\n",
    "    print(f\"  watertight:             {watertight}\")\n",
    "    print(f\"  orientable:             {orientable}\")\n",
    "\n",
    "    geoms = [mesh]\n",
    "    if not edge_manifold:\n",
    "        edges = mesh.get_non_manifold_edges(allow_boundary_edges=True)\n",
    "        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 0)))\n",
    "    if not edge_manifold_boundary:\n",
    "        edges = mesh.get_non_manifold_edges(allow_boundary_edges=False)\n",
    "        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (0, 1, 0)))\n",
    "    if not vertex_manifold:\n",
    "        verts = np.asarray(mesh.get_non_manifold_vertices())\n",
    "        pcl = o3d.geometry.PointCloud(\n",
    "            points=o3d.utility.Vector3dVector(np.asarray(mesh.vertices)[verts]))\n",
    "        pcl.paint_uniform_color((0, 0, 1))\n",
    "        geoms.append(pcl)\n",
    "    if self_intersecting:\n",
    "        intersecting_triangles = np.asarray(\n",
    "            mesh.get_self_intersecting_triangles())\n",
    "        intersecting_triangles = intersecting_triangles[0:1]\n",
    "        intersecting_triangles = np.unique(intersecting_triangles)\n",
    "        print(\"  # visualize self-intersecting triangles\")\n",
    "        triangles = np.asarray(mesh.triangles)[intersecting_triangles]\n",
    "        edges = [\n",
    "            np.vstack((triangles[:, i], triangles[:, j]))\n",
    "            for i, j in [(0, 1), (1, 2), (2, 0)]\n",
    "        ]\n",
    "        edges = np.hstack(edges).T\n",
    "        edges = o3d.utility.Vector2iVector(edges)\n",
    "        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 1)))\n",
    "    o3d.visualization.draw_geometries(geoms, mesh_show_back_face=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check_properties('Knot', o3dtut.get_knot_mesh())\n",
    "check_properties('Moebius', o3d.geometry.TriangleMesh.create_moebius(twists=1))\n",
    "check_properties(\"non-manifold edge\", o3dtut.get_non_manifold_edge_mesh())\n",
    "check_properties(\"non-manifold vertex\", o3dtut.get_non_manifold_vertex_mesh())\n",
    "check_properties(\"open box\", o3dtut.get_open_box_mesh())\n",
    "check_properties(\"intersecting_boxes\", o3dtut.get_intersecting_boxes_mesh())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh filtering\n",
    "Open3D contains a number of methods to filter meshes. In the following we show the implemented filters to smooth noisy triangle meshes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average filter\n",
    "The simplest filter is the average filter. A given vertex $v_i$ is given by the average of the adjacent vertices $\\mathcal{N}$\n",
    "\n",
    "\\begin{equation}\n",
    "v_i = \\frac{v_i + \\sum_{n \\in \\mathcal{N}} v_n}{|N| + 1} \\,.\n",
    "\\end{equation}\n",
    "\n",
    "This filter can be used to denoise meshes as demonstrated in the code below. The parameter `number_of_iterations` in the function `filter_smooth_simple` defines the how often the filter is applied to the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('create noisy mesh')\n",
    "mesh_in = o3dtut.get_knot_mesh()\n",
    "vertices = np.asarray(mesh_in.vertices)\n",
    "noise = 5\n",
    "vertices += np.random.uniform(0, noise, size=vertices.shape)\n",
    "mesh_in.vertices = o3d.utility.Vector3dVector(vertices)\n",
    "mesh_in.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_in])\n",
    "\n",
    "print('filter with average with 1 iteration')    \n",
    "mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=1)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])\n",
    "\n",
    "print('filter with average with 5 iterations')    \n",
    "mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=5)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplacian\n",
    "Another important mesh filter is the Laplacian defined as \n",
    "\n",
    "\\begin{equation}\n",
    "v_i = v_i \\cdot \\lambda \\sum_{n \\in N} w_n v_n - v_i \\,,\n",
    "\\end{equation}\n",
    "\n",
    "where $\\lambda$ is the strenght of the filter and $w_n$ are normalized weights that relate to the distance of the neighbouring vertices. The filter is implemented in `filter_smooth_laplacian` and has the parameters `number_of_iterations` and `lambda`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('filter with Laplacian with 10 iterations')    \n",
    "mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=10)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])\n",
    "\n",
    "print('filter with Laplacian with 50 iterations')    \n",
    "mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=50)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Taubin filter\n",
    "The problem with the average and Laplacian filter is that they lead to a shrinkage of the triangle mesh. [\\[Taubin1995\\]](../reference.html#Taubin1995) showed that the application of two Laplacian filters with different $\\lambda$ parameters can prevent the mesh shrinkage. The filter is implemented in `filter_smooth_taubin`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('filter with Taubin with 10 iterations')  \n",
    "mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=10)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])\n",
    "\n",
    "print('filter with Taubin with 100 iterations')  \n",
    "mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=100)\n",
    "mesh_out.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh_out])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling\n",
    "Open3D includes functions to sample point clouds from a triangle mesh. The simplest method is `sample_points_uniformly` that uniformly samples points from the 3D surface based on the triangle area. The parameter `number_of_points` defines how many points are sampled from the triangle surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3d.geometry.TriangleMesh.create_sphere()\n",
    "mesh.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh])\n",
    "pcd = mesh.sample_points_uniformly(number_of_points=500)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3dtut.get_bunny_mesh()\n",
    "mesh.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh])\n",
    "pcd = mesh.sample_points_uniformly(number_of_points=500)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uniform sampling can yield clusters of points on the surface. Poisson disk sampling can evenly distributes the points on the surface. The method `sample_points_poisson_disk` implements sample elimination. Therefore, the method starts with a sampled point cloud and removes point to satisfy the sampling criterion. The method supports two options to provide the initial point cloud')\n",
    "\n",
    "1. Default via the parameter `init_factor`: The method first samples uniformly a point cloud from the mesh with `init_factor` x `number_of_points` and uses this for the elimination\n",
    "2. One can provide a point cloud and pass it to the `sample_points_poisson_disk` method. Then, this point cloud is used for elimination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3d.geometry.TriangleMesh.create_sphere()\n",
    "pcd = mesh.sample_points_poisson_disk(number_of_points=500, init_factor=5)\n",
    "o3d.visualization.draw_geometries([pcd])\n",
    "\n",
    "pcd = mesh.sample_points_uniformly(number_of_points=2500)\n",
    "pcd = mesh.sample_points_poisson_disk(number_of_points=500, pcl=pcd)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3dtut.get_bunny_mesh()\n",
    "pcd = mesh.sample_points_poisson_disk(number_of_points=500, init_factor=5)\n",
    "o3d.visualization.draw_geometries([pcd])\n",
    "\n",
    "pcd = mesh.sample_points_uniformly(number_of_points=2500)\n",
    "pcd = mesh.sample_points_poisson_disk(number_of_points=500, pcl=pcd)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh subdivision\n",
    "In mesh subdivision we divide each triangles into a number of smaller triangles. In the simplest case, we compute the midpoint of each side per triangle and divide the triangle into four smaller triangles. This is implemented in the `subdivide_midpoint` function. The 3D surface and area stays the same, but the number of vertices and triangles increases. The parameter `number_of_iterations` defines how often this process is repeated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3d.geometry.TriangleMesh.create_box()\n",
    "mesh.compute_vertex_normals()\n",
    "print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)\n",
    "mesh = mesh.subdivide_midpoint(number_of_iterations=1)\n",
    "print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open3D implements an additional subdivision method based on [\\[Loop1987\\]](../reference.html#Loop1987). The method is based on a quartic box spline, which generate $C^2$ continuous limit surfaces everywhere except at extraordinary vertices where they are $C^1$ continuous. This leads to smoother corners."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3d.geometry.TriangleMesh.create_sphere()\n",
    "mesh.compute_vertex_normals()\n",
    "print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)\n",
    "mesh = mesh.subdivide_loop(number_of_iterations=2)\n",
    "print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3dtut.get_knot_mesh()\n",
    "mesh.compute_vertex_normals()\n",
    "print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)\n",
    "mesh = mesh.subdivide_loop(number_of_iterations=1)\n",
    "print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh simplification\n",
    "Sometimes we want to represent a high-resolution mesh with fewer triangles and vertices, but the low-resolution mesh should still be close to the high-resolution mesh. For this purpose Open3D implements a number of mesh simplification methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vertex clustering\n",
    "The vertex clustering method pools all vertices that fall into a voxel of a given size to a single vertex. The method is implemented in `simplify_vertex_clustering` and has as parameters `voxel_size` that defines the size of the voxel grid and `contraction` that defines how the vertices are pooled. `o3d.geometry.SimplificationContraction.Average` computes a simple average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_in = o3dtut.get_bunny_mesh()\n",
    "print(f'Input mesh has {len(mesh_in.vertices)} vertices and {len(mesh_in.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh_in])\n",
    "\n",
    "voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 32\n",
    "print(f'voxel_size = {voxel_size:e}')\n",
    "mesh_smp = mesh_in.simplify_vertex_clustering(\n",
    "    voxel_size=voxel_size,\n",
    "    contraction=o3d.geometry.SimplificationContraction.Average)\n",
    "print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh_smp])\n",
    "\n",
    "voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 16\n",
    "print(f'voxel_size = {voxel_size:e}')\n",
    "mesh_smp = mesh_in.simplify_vertex_clustering(\n",
    "    voxel_size=voxel_size,\n",
    "    contraction=o3d.geometry.SimplificationContraction.Average)\n",
    "print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh_smp])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh decimation\n",
    "Another catogory of mesh simplification methods is mesh decimation that operates in incremental steps. We select a single triangle that minimizes a error metric and remove it. This is repeated until a required number of triangles is achieved. Open3D implements `simplify_quadric_decimation` that minimizes error quadrics (distances to neighbouring planes). The parameter `target_number_of_triangles` defines the stopping critera of the decimation algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_smp = mesh_in.simplify_quadric_decimation(\n",
    "    target_number_of_triangles=6500)\n",
    "print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh_smp])\n",
    "\n",
    "mesh_smp = mesh_in.simplify_quadric_decimation(\n",
    "    target_number_of_triangles=1700)\n",
    "print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')\n",
    "o3d.visualization.draw_geometries([mesh_smp])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connected components\n",
    "The result of various reconstruction methods (e.g., [RGBD Integration](../Advanced/rgbd_integration.ipynb) is not always a single triangle mesh, but a number of meshes. Some of the smaller parts are due to noise and we most likely want to remove them. Open3D implements a connected components algorithm `cluster_connected_triangles` that assigns each triangle to a cluster of connected triangles. It returns for each triangle the index of the cluster in `triangle_clusters`, and per cluster the number of triangles in `cluster_n_triangles` and the surface area of the cluster in `cluster_area`.\n",
    "\n",
    "The code below shows the application of `cluster_connected_triangles` and how it can be used to remove spurious triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Generate data\")\n",
    "mesh = o3dtut.get_bunny_mesh().subdivide_midpoint(number_of_iterations=2)\n",
    "vert = np.asarray(mesh.vertices)\n",
    "min_vert, max_vert = vert.min(axis=0), vert.max(axis=0)\n",
    "for _ in range(30):\n",
    "    cube = o3d.geometry.TriangleMesh.create_box()\n",
    "    cube.scale(0.005, center=cube.get_center())\n",
    "    cube.translate(\n",
    "        (\n",
    "            np.random.uniform(min_vert[0], max_vert[0]),\n",
    "            np.random.uniform(min_vert[1], max_vert[1]),\n",
    "            np.random.uniform(min_vert[2], max_vert[2]),\n",
    "        ),\n",
    "        relative=False,\n",
    "    )\n",
    "    mesh += cube\n",
    "mesh.compute_vertex_normals()\n",
    "print(\"Show input mesh\")\n",
    "o3d.visualization.draw_geometries([mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Cluster connected triangles\")\n",
    "with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug) as cm:\n",
    "    triangle_clusters, cluster_n_triangles, cluster_area = (\n",
    "        mesh.cluster_connected_triangles())\n",
    "triangle_clusters = np.asarray(triangle_clusters)\n",
    "cluster_n_triangles = np.asarray(cluster_n_triangles)\n",
    "cluster_area = np.asarray(cluster_area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Show mesh with small clusters removed\")\n",
    "mesh_0 = copy.deepcopy(mesh)\n",
    "triangles_to_remove = cluster_n_triangles[triangle_clusters] < 100\n",
    "mesh_0.remove_triangles_by_mask(triangles_to_remove)\n",
    "o3d.visualization.draw_geometries([mesh_0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Show largest cluster\")\n",
    "mesh_1 = copy.deepcopy(mesh)\n",
    "largest_cluster_idx = cluster_n_triangles.argmax()\n",
    "triangles_to_remove = triangle_clusters != largest_cluster_idx\n",
    "mesh_1.remove_triangles_by_mask(triangles_to_remove)\n",
    "o3d.visualization.draw_geometries([mesh_1])"
   ]
  }
 ],
 "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": 2
}