docs/gen_3d_mesh.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "painful-march",
   "metadata": {},
   "source": [
    "# Generate 3D mesh from point cloud\n",
    "\n",
    "Based on the article [_\"5-Step Guide to generate 3D meshes from point\n",
    "clouds with Python\"_\n",
    "](https://towardsdatascience.com/5-step-guide-to-generate-3d-meshes-from-point-clouds-with-python-36bad397d8ba) by Florent Poux. In this notebook we will generate a 3D mesh for the Kermit dataset processed with [`linux-photogrammetry-tools`](https://github.com/epassaro/linux-photogrammetry-tools), so some steps could differ from the original source.\n",
    "\n",
    "### Prerequisites\n",
    "\n",
    "To run this notebook create a new `conda` environment:\n",
    "\n",
    "```\n",
    "$ conda create -n open3d -c open3d-admin -c conda-forge open3d=0.12\n",
    "```\n",
    "\n",
    "## 1. Load point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surgical-reduction",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import open3d as o3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "thick-circle",
   "metadata": {},
   "outputs": [],
   "source": [
    "point_cloud = np.loadtxt(os.path.join(\"../examples/kermit/output/option-0000.ply\"), skiprows=14)\n",
    "output_path = os.path.join(\".\", \"output\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "urban-outline",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcd = o3d.geometry.PointCloud()\n",
    "pcd.points = o3d.utility.Vector3dVector(point_cloud[:,:3])\n",
    "pcd.normals = o3d.utility.Vector3dVector(point_cloud[:,3:6])\n",
    "pcd.colors = o3d.utility.Vector3dVector(point_cloud[:,6:-1]/255)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "removable-african",
   "metadata": {},
   "outputs": [],
   "source": [
    "#o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recovered-latin",
   "metadata": {},
   "source": [
    "## 2.1 Meshing with ball pivoting algorithm (BPA)\n",
    "\n",
    "We first compute the necessary radius parameter based on the average distances computed from all the distances between points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "proprietary-expression",
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = pcd.compute_nearest_neighbor_distance()\n",
    "avg_dist = np.mean(distances)\n",
    "radius = 3*avg_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "friendly-allergy",
   "metadata": {},
   "source": [
    "In one command line, we can then create a mesh and store it in the `bpa_mesh` variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "critical-selling",
   "metadata": {},
   "outputs": [],
   "source": [
    "bpa_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(pcd,o3d.utility.DoubleVector([radius, radius * 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spoken-vacation",
   "metadata": {},
   "source": [
    "Before exporting the mesh, we can downsample the result to an acceptable number of triangles, for example, 100k triangles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "swedish-document",
   "metadata": {},
   "outputs": [],
   "source": [
    "dec_mesh = bpa_mesh.simplify_quadric_decimation(100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "working-october",
   "metadata": {},
   "source": [
    "Additionally, if you think the mesh can present some weird artifacts, you can run the following commands to ensure its consistency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unlike-ottawa",
   "metadata": {},
   "outputs": [],
   "source": [
    "dec_mesh.remove_degenerate_triangles()\n",
    "dec_mesh.remove_duplicated_triangles()\n",
    "dec_mesh.remove_duplicated_vertices()\n",
    "dec_mesh.remove_non_manifold_edges()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "functional-whale",
   "metadata": {},
   "source": [
    "## 2.2 Meshing with Poisson surface reconstruction\n",
    "\n",
    "To get results with Poisson, it is very straightforward. You just have to adjust the\n",
    "parameters that you pass to the function as described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "close-offset",
   "metadata": {},
   "outputs": [],
   "source": [
    "poisson_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, depth=8, width=0, scale=1.1, \\\n",
    "                                                                         linear_fit=False)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "detected-negative",
   "metadata": {},
   "source": [
    "To get a clean result, it is often necessary to add a cropping step to clean unwanted artifacts. For this, we compute the initial bounding-box containing the raw point cloud, and we\n",
    "use it to filter all surfaces from the mesh outside the bounding-box:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "exterior-jones",
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = pcd.get_axis_aligned_bounding_box()\n",
    "p_mesh_crop = poisson_mesh.crop(bbox)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "violent-spirit",
   "metadata": {},
   "source": [
    "## 3. Output and visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fiscal-greenhouse",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(output_path, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sunset-technical",
   "metadata": {},
   "outputs": [],
   "source": [
    "o3d.io.write_triangle_mesh(os.path.join(output_path, \"bpa_mesh.ply\"), dec_mesh)\n",
    "o3d.io.write_triangle_mesh(os.path.join(output_path, \"p_mesh.ply\"), p_mesh_crop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "decreased-contract",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lod_mesh_export(mesh, lods, extension, path):\n",
    "    \n",
    "    mesh_lods={}\n",
    "    for i in lods:\n",
    "        mesh_lod = mesh.simplify_quadric_decimation(i)\n",
    "        o3d.io.write_triangle_mesh(os.path.join(output_path, f\"lod_{str(i)}{extension}\"), mesh_lod)\n",
    "        mesh_lods[i] = mesh_lod\n",
    "    \n",
    "    print(\"generation of \"+str(i)+\" LoD successful\")\n",
    "    \n",
    "    return mesh_lods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fixed-honduras",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_lods = lod_mesh_export(p_mesh_crop, [100000, 50000, 10000, 1000, 100], \".ply\", output_path)\n",
    "#my_lods2 = lod_mesh_export(bpa_mesh, [8000, 800, 300], \".ply\", output_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "wrong-salmon",
   "metadata": {},
   "outputs": [],
   "source": [
    "o3d.visualization.draw_geometries([my_lods[100000]])"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}