{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Работа с полигональными моделями и облаками точек. MeshLab / `pymeshlab`\n",
    "\n",
    "Полигональная/фасетная модель (меш) это набор многоугольников, обычно треугольников. Задается в виде списка точек (которые задаются в виде списка координат) и еще троек.\n",
    "\n",
    "На данные, что на входные, что на выходные, очень полезно смотреть глазами.\n",
    "\n",
    "С картинками, текстами и аудио это делается тривиально.\n",
    "\n",
    "Рассмотрим, как это делается с 3D данными.\n",
    "\n",
    "\n",
    "[Download MeshLab GUI](https://www.meshlab.net/#download)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "< Секция с демонстрацией открытия меша, семплингом, подсчетом метрик >"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install open3d\n",
    "# !pip install pymeshlab\n",
    "# !pip install --upgrade PyMCubes\n",
    "# !pip install edt\n",
    "# !pip install mesh2sdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Что происходит внутри"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 320\n",
    "h = 200"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table align='left'>\n",
    "<tr>\n",
    "<td><img src='imgs/mesh.png' width=str(w) height=str(h)/></td>\n",
    "<td><img src='imgs/pc.png' width=str(w) height=str(h)/></td>\n",
    "<td><img src='imgs/m.png' width=str(w) height=str(h)/></td>\n",
    "</tr>\n",
    "\n",
    "<table align='left'>\n",
    "<tr>\n",
    "<td><img src='imgs/metric.png' width=str(w) height=str(h)/></td>\n",
    "<td><img src='imgs/geo.png' width=str(w) height=str(h)/></td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Как то же самое сделать программно"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import open3d as o3d\n",
    "import pymeshlab as pml\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"Легкая\" Визуализация удобнее в Open3D\n",
    "\n",
    "armadillo_mesh = o3d.data.ArmadilloMesh() # скачиваем из интернета мешфайл\n",
    "path = armadillo_mesh.path\n",
    "\n",
    "# path = 'armadillo.obj'\n",
    "mesh = o3d.io.read_triangle_mesh(path)\n",
    "mesh.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Вычисления/Операции над мешами делаются в MeshLab\n",
    "\n",
    "ms = pml.MeshSet() # создаем мешсет\n",
    "ms.load_new_mesh(path)\n",
    "#https://pymeshlab.readthedocs.io/en/latest/filter_list.html#generate_sampling_poisson_disk\n",
    "ms.generate_sampling_poisson_disk(samplenum = 10000) # создается объект ms[1], в котором только вершины\n",
    "ms.generate_surface_reconstruction_screened_poisson() # создается объект ms[2], это меш\n",
    "ms.compute_scalar_by_distance_from_another_mesh_per_vertex(measuremesh=0, refmesh=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scalar_np = ms[0].vertex_scalar_array()\n",
    "print('per-vertex scalar values are:', scalar_np)\n",
    "vertices_np = ms[0].vertex_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize_array(arr:np.array):\n",
    "    '''     remaps any array to [0,1]    '''\n",
    "    L = arr.max() - arr.min()\n",
    "    return (arr - arr.min()) / L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcd = o3d.geometry.PointCloud()\n",
    "pcd.points = o3d.utility.Vector3dVector(vertices_np)\n",
    "normalized = normalize_array(scalar_np)\n",
    "cmap = plt.cm.RdBu\n",
    "colors = cmap(normalized)[:,:3]\n",
    "pcd.colors = o3d.utility.Vector3dVector(colors)\n",
    "\n",
    "# открывается отдельное окно\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Signed Distance Function\n",
    "Определение в стиле Computer Science / Computer Graphics\n",
    "\n",
    "\\begin{equation}\n",
    "    SDF(x) = \\left\\{\\begin{array}{rl}\n",
    "        dist(x, \\partial\\Omega), &\\textrm{если $x$ снаружи объекта } \\Omega, \\\\\n",
    "        -dist(x, \\partial\\Omega), &\\textrm{иначе,} \\\\\n",
    "    \\end{array}\\right.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Определение в стиле математики\n",
    "\\begin{equation}\n",
    "\\begin{array}{rl}\n",
    "    | \\nabla u(x)| &= 1, x \\in \\mathbb{R}^3 \\\\\n",
    "    u|_{\\partial \\Omega} &= 0 \\\\\n",
    "    u(x) >0 &\\textrm{если $x$ снаружи объекта }\n",
    "\\end{array}\\\n",
    "\\end{equation}\n",
    "\n",
    "(Сайт для доп. изучения)(https://iquilezles.org/articles/distfunctions/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Также важное понятие:\n",
    "\n",
    "\n",
    "Усечённая знаковая функция расстояния (Truncated Signed Distance Function)\n",
    "\\begin{equation}\n",
    "    TSDF(p) = \\left\\{\\begin{array}{rl}\n",
    "        \\mu, &\\textrm{if } \\mu \\leq SDF(p) \\\\\n",
    "        SDF(p), &\\textrm{if } -\\mu < SDF(p) < \\mu \\\\\n",
    "        -\\mu, &\\textrm{if } SDF(p) \\leq -\\mu. \\\\\n",
    "    \\end{array}\\right.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\n",
    "1. Исторически важные статьи TSDF Fusion и KinectFusion\n",
    "SOTA донейронной эпохи. Работает в реальном времени на GPU времен 2010г.\n",
    "Идея была такая, что по карте глубины можно по формуле считать \"проективный SDF\", который (примерно) верен вблизи поверхности. Если поснимать объект сразных сторон и усреднить, вблизи поверхности получается очень близко к настоящему SDF. Но нужно обрезать результат в окрестности от поверхности, чтобы не начать ловить артефакты.\n",
    "\n",
    "2. В современной литературе даже когда слова Truncated нет в названии, в статье де-факто используется именно TSDF.\n",
    "\n",
    "\"Честный\" SDF можно получить *только* \n",
    "\n",
    "- для фигур, для которых он известен аналитически (в том числе для меш, но это затратно)\n",
    "- честным решателем уравнения эйконала, таким как ```edt``` (есть методы быстрее чем N^3)\n",
    "- через фитинг эйконал-лосса на нейронную сеть или еще куда\n",
    "\n",
    "$$L = (\\left| \\nabla u_\\theta(x) \\right| -1 )^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D SDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table align='left'>\n",
    "    <tr>\n",
    "    <img src='imgs/sdf2d.jpg' width=440 height=str(h)/>\n",
    "    <img src='imgs/sdf4.png' width=440 height=str(h)/>\n",
    "    </tr>\n",
    "    <table align='left'>\n",
    "</table>\n",
    "\n",
    "<table align='left'>\n",
    "    <tr>\n",
    "    <td><img src='imgs/u1.png' width=str(w) height=str(h)/></td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D SDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ModelNet plane-0629 (watertight version), SDF, TSDF (normalized)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table align='left'>\n",
    "    <tr>\n",
    "    <img src='imgs/dat_plane00s.png' width=440 height=str(h)/>\n",
    "    </tr>\n",
    "    <tr>\n",
    "    <img src='imgs/sdf_full.png' width=str(w) height=str(h)/>\n",
    "    </tr>\n",
    "    <tr>\n",
    "    <img src='imgs/tsdf_full.png' width=str(w) height=str(h)/>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mcubes\n",
    "import edt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128\n",
    "x, y, z = np.mgrid[:N, :N, :N]\n",
    "\n",
    "binary_sphere = (x - N//2)**2 + (y - N//2)**2 + (z - N//2)**2 - (N//4)**2 < 0\n",
    "vertices, triangles = mcubes.marching_cubes(binary_sphere, 0.5) # 1 - binary_sphere to keep normals outwards\n",
    "mcubes.export_obj(vertices, triangles, \"sphere.obj\", \"MySphere\")\n",
    "\n",
    "smoothed_sphere = mcubes.smooth(binary_sphere)\n",
    "vertices_smooth, triangles_smooth = mcubes.marching_cubes(smoothed_sphere, 0)\n",
    "mcubes.export_obj(vertices_smooth, triangles_smooth, \"sphere_smooth.obj\", \"MySphere\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(binary_sphere[N//2,:,:])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sdf(x):\n",
    "    return (1-x)*(edt.edt(1-x)-1) - edt.edt(x)*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf_data = sdf(binary_sphere)\n",
    "vertices_sdf, triangles_sdf = mcubes.marching_cubes(sdf_data, 0)\n",
    "mcubes.export_obj(vertices_sdf, triangles_sdf, \"sphere_from_sdf.obj\")\n",
    "\n",
    "data = sdf_data[N//2, :, :]\n",
    "plt.imshow(data)\n",
    "plt.colorbar()\n",
    "plt.contour(data, levels=[0], colors=['white'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = o3d.io.read_triangle_mesh('chair.obj')\n",
    "mesh.compute_vertex_normals()\n",
    "o3d.visualization.draw_geometries([mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [Open3D INFO]   -- Render mode control --\n",
    "# [Open3D INFO]     L            : Turn on/off lighting.\n",
    "# [Open3D INFO]     +/-          : Increase/decrease point size.\n",
    "# [Open3D INFO]     Ctrl + +/-   : Increase/decrease width of geometry::LineSet.\n",
    "# [Open3D INFO]     N            : Turn on/off point cloud normal rendering.\n",
    "# [Open3D INFO]     S            : Toggle between mesh flat shading and smooth shading.\n",
    "# [Open3D INFO]     W            : Turn on/off mesh wireframe.\n",
    "# [Open3D INFO]     B            : Turn on/off back face rendering.\n",
    "# [Open3D INFO]     I            : Turn on/off image zoom in interpolation.\n",
    "# [Open3D INFO]     T            : Toggle among image render:\n",
    "# [Open3D INFO]                    no stretch / keep ratio / freely stretch.\n",
    "# [Open3D INFO]     Ctrl + 0..4,9: Set mesh color option.\n",
    "# [Open3D INFO]                    0 - Default behavior, render uniform gray color.\n",
    "# [Open3D INFO]                    1 - Render point color.\n",
    "# [Open3D INFO]                    2 - x coordinate as color.\n",
    "# [Open3D INFO]                    3 - y coordinate as color.\n",
    "# [Open3D INFO]                    4 - z coordinate as color.\n",
    "# [Open3D INFO]                    9 - normal as color.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "# import trimesh\n",
    "import mesh2sdf\n",
    "\n",
    "import time\n",
    "\n",
    "filename = sys.argv[1] if len(sys.argv) > 1 else  \\\n",
    "    os.path.join(os.path.dirname(__file__), 'data', 'plane.obj')\n",
    "\n",
    "mesh_scale = 0.8\n",
    "size = 128\n",
    "level = 2 / size\n",
    "\n",
    "# mesh = trimesh.load('chair.obj')\n",
    "\n",
    "# normalize mesh\n",
    "vertices = np.asarray(mesh.vertices)\n",
    "faces = np.asarray(mesh.triangles)\n",
    "bbmin = vertices.min(0)\n",
    "bbmax = vertices.max(0)\n",
    "center = (bbmin + bbmax) * 0.5\n",
    "scale = 2.0 * mesh_scale / (bbmax - bbmin).max()\n",
    "vertices = (vertices - center) * scale\n",
    "\n",
    "# fix mesh\n",
    "t0 = time.time()\n",
    "sdf, mesh = mesh2sdf.compute(\n",
    "    vertices, faces, size, fix=True, level=level, return_mesh=True)\n",
    "t1 = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sdf[size//2, :, :]\n",
    "plt.imshow(data)\n",
    "plt.contour(data, levels=[0], colors=['white'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Конверсия в бинарное:\n",
    "Сам по себе набор треугольников не несет информации о том, лежит ли та или иная точка внутри или снаружи.\n",
    "\n",
    "Старая техника получения сдф\n",
    "1) из меша можно получить TSDF, рендеря карту глубины с разных камер.\n",
    "2) TSDF будет очень хорошо аппроксимировать в окрестности поверхности\n",
    "3) каким-нибудь точным солвером уравнения эйконала (воксельным) найти SDF, склеить их вместе"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Курс \"Трехмерное компьютерное зрение\" в ШАД, А.Конушин и команда.\n",
    "Семинар подготовил А.Бойко."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "shad2",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}