mesh_creation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Warning] Since Open3D 0.15, installing Open3D via conda is deprecated. Please re-install Open3D via: `pip install open3d -U`.\n",
      "Jupyter environment detected. Enabling Open3D WebVisualizer.\n",
      "[Open3D INFO] WebRTC GUI backend enabled.\n",
      "[Open3D INFO] WebRTCWindowSystem: HTTP handshake server disabled.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pyvista as pv\n",
    "import open3d\n",
    "import modeling\n",
    "from modeling import *\n",
    "import scipy.ndimage\n",
    "#import scipy.stats\n",
    "import copy\n",
    "from utils import *\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import cm\n",
    "import tensorflow as tf\n",
    "\n",
    "\n",
    "\n",
    "silver = [0.753, 0.753, 0.753]#{'color': [0.753, 0.753, 0.753], 'material': \n",
    "a_si = [0.600, 0.460, 0.357]\n",
    "gold = [212/255, 175/255, 55/255]\n",
    "\n",
    "si_material = open3d.visualization.Material('defaultLitSSR')\n",
    "si_material.scalar_properties['metallic'] = 0.0\n",
    "si_material.scalar_properties['roughness'] = 0.218\n",
    "si_material.vector_properties['base_color'] = [*a_si, 1]\n",
    "ag_material = open3d.visualization.Material('defaultLitSSR')\n",
    "ag_material.scalar_properties['metallic'] = 1.0\n",
    "ag_material.scalar_properties['roughness'] = 0.666\n",
    "ag_material.vector_properties['base_color'] = [*silver, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'L': 256, 'theta': 85.5, 'phi': 0, 'H': 64, 'D': 10, 'turns': 0.85, 'species': [1], 'weights': array([[ 1.        ,  0.        ,  0.        ],\n",
      "       [ 0.        ,  1.        ,  0.        ],\n",
      "       [ 0.2       , -0.30000001,  3.        ]]), 'repetition': 524288, 'time': 3105.510812997818}, {'L': 256, 'theta': 85.5, 'phi': 306.0, 'H': 96, 'D': 30, 'turns': 0.15, 'species': [2], 'weights': array([[ 1.        ,  0.        ,  0.        ],\n",
      "       [ 0.        ,  1.        ,  0.        ],\n",
      "       [ 0.2       , -0.30000001,  3.        ]]), 'repetition': 65536, 'time': 659.2555365562439}]\n"
     ]
    }
   ],
   "source": [
    "filename = 'SCTF_Si_Ag_L256_x1.0_Th85.5_D30_N589824_1659041585'\n",
    "\n",
    "grid, deposited, params = loadSparse('structures//' + filename + '.npz')\n",
    "print(params)\n",
    "\n",
    "grid_full = matrixFromSparse(grid)\n",
    "grid_si = np.where(grid_full==1, 1, 0)\n",
    "grid_ag = np.where(grid_full==2, 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open3D works better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'L': 256, 'theta': 85.5, 'phi': 0, 'H': 64, 'D': 10, 'turns': 0.85, 'species': [1], 'weights': array([[ 1.        ,  0.        ,  0.        ],\n",
      "       [ 0.        ,  1.        ,  0.        ],\n",
      "       [ 0.2       , -0.30000001,  3.        ]]), 'repetition': 524288, 'time': 3105.510812997818}, {'L': 256, 'theta': 85.5, 'phi': 306.0, 'H': 96, 'D': 30, 'turns': 0.15, 'species': [2], 'weights': array([[ 1.        ,  0.        ,  0.        ],\n",
      "       [ 0.        ,  1.        ,  0.        ],\n",
      "       [ 0.2       , -0.30000001,  3.        ]]), 'repetition': 65536, 'time': 659.2555365562439}]\n",
      "Disconnected floating components.\n",
      "452 Ag regions:\n",
      "\t75 able to mesh; 377 unable for 4122 points left out.\n",
      "190 Si regions:\n",
      "\t94 able to mesh; 96 unable for 1052 points left out.\n",
      "Created surface meshes.\n"
     ]
    }
   ],
   "source": [
    "filename = 'STF_Si_Ag_L768_x1.0_Th85.5_D60_N9437184_1660238350'\n",
    "filename = 'SCTF_Si_Ag_L256_x1.0_Th85.5_D30_N589824_1659041585'\n",
    "si_meshes, ag_meshes = reconstructSiAgMeshColumns(filename, start=6, color_by_column=True, remove_floating=True)\n",
    "open3d.visualization.draw_geometries([*ag_meshes, *si_meshes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Disconnected floating components.\n",
      "Hollowed components.\n",
      "Created point clouds.\n",
      "Created surface meshes.\n"
     ]
    }
   ],
   "source": [
    "''' Clouds all at once '''\n",
    "filename = 'STF_Si_Ag_L768_x1.0_Th85.5_D60_N9437184_1660238350'\n",
    "si_mesh, ag_mesh = modeling.reconstructSiAgMeshesFull('SCTF_Si_Ag_L256_x1.0_Th85.5_D30_N589824_1659041585')\n",
    "open3d.visualization.draw_geometries([ag_mesh, si_mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Out of 125:\n",
      "\t43 not orientable\n",
      "\t4 are not vertex_manifold\n",
      "\t43 are not edge manifold\n",
      "\t0 are self intersecting\n",
      "\t0 else\n"
     ]
    }
   ],
   "source": [
    "#open3d.visualization.draw_geometries([si_meshes[-1]])\n",
    "#open3d.visualization.draw_geometries([*si_meshes, *ag_meshes, *si_comps, *ag_comps])\n",
    "nor_list = []\n",
    "van_list = []\n",
    "ean_list = []\n",
    "sin_list = []\n",
    "els_list = []\n",
    "for i in range(len(ag_meshes)):\n",
    "    ag_meshes[i].remove_degenerate_triangles()\n",
    "    if not ag_meshes[i].is_orientable():\n",
    "        ag_meshes[i].paint_uniform_color([1, 0.5, 0.5])\n",
    "        nor_list.append(i)\n",
    "    if not ag_meshes[i].is_vertex_manifold():\n",
    "        ag_meshes[i].paint_uniform_color([0.5, 1, 1]) # teal\n",
    "        van_list.append(i)\n",
    "    if not ag_meshes[i].is_edge_manifold():\n",
    "        ag_meshes[i].paint_uniform_color([1, 0.5, 1]) # purple\n",
    "        ean_list.append(i)\n",
    "    if ag_meshes[i].is_self_intersecting():\n",
    "        ag_meshes[i].paint_uniform_color([0, 1, 0.5]) # \n",
    "        sin_list.append(i)\n",
    "    #else:\n",
    "    #    ag_meshes[i].paint_uniform_color(silver)\n",
    "    #    els_list.append(i)\n",
    "print(\"Out of {}:\\n\\t{} not orientable\\n\\t{} are not vertex_manifold\\n\\t{} are not edge manifold\\n\\t{} are self intersecting\\n\\t{} else\".\n",
    "      format(len(ag_meshes), len(nor_list), len(van_list), len(ean_list), len(sin_list), len(els_list)))\n",
    "open3d.visualization.draw_geometries([*si_meshes[:-1], *ag_meshes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[180 182]]\n",
      "19\n"
     ]
    }
   ],
   "source": [
    "index = 1\n",
    "lst = van_list\n",
    "mesh = ag_meshes[lst[index]]\n",
    "edges = mesh.get_non_manifold_edges(allow_boundary_edges=False)\n",
    "print(np.asarray(edges))\n",
    "line = open3d.geometry.LineSet(points=mesh.vertices, lines=edges)\n",
    "line.colors = open3d.utility.Vector3dVector([[0,1,0]])\n",
    "mesh.vertex_colors[edges[0][0]] = np.array([0, 1, 0])\n",
    "mesh.vertex_colors[edges[0][1]] = np.array([0, 1, 0])\n",
    "print(lst[index])\n",
    "open3d.visualization.draw_geometries([ag_comps[lst[index]], line])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "render = open3d.visualization.rendering.OffscreenRenderer(1024, 1024)\n",
    "silvermat = open3d.visualization.rendering.MaterialRecord()\n",
    "silvermat.base_color = [*gold, 1]\n",
    "silvermat.shader = 'defaultLit'\n",
    "silvermat.base_metallic = 2\n",
    "siliconmat = open3d.visualization.rendering.MaterialRecord()\n",
    "siliconmat.base_color = [*a_si, 1]\n",
    "siliconmat.shader = 'defaultLit'\n",
    "siliconmat.base_roughness = 16\n",
    "\n",
    "i = 0\n",
    "for m in si_meshes:\n",
    "    if m != si_meshes[-1]:\n",
    "        render.scene.add_geometry(\"Si_{}\".format(i), m, siliconmat)\n",
    "        i += 1\n",
    "i = 0\n",
    "for m in ag_meshes:\n",
    "    render.scene.add_geometry(\"Ag_{}\".format(i), m, silvermat)\n",
    "    i += 1\n",
    "render.setup_camera(40, [384, 384, 0], [-384, 384, 512], [0, 0, 1])\n",
    "render.scene.scene.set_sun_light([0, -0.5, -1], [1, 1, 1], 7500)\n",
    "render.scene.scene.enable_sun_light(True)\n",
    "render.scene.show_axes(False)\n",
    "#render.set_clear_color([0,0,0,1])\n",
    "img = render.render_to_image()\n",
    "open3d.io.write_image(\"models/\" + filename + \"_front.png\", img)\n",
    "render.setup_camera(40, [384, 384, 0], [0, 0, 512], [0, 0, 1])\n",
    "img = render.render_to_image()\n",
    "open3d.io.write_image(\"models/\" + filename + \"_diagonal.png\", img)\n",
    "render.setup_camera(40, [384, 384, 0], [384, 384, 768], [0, 0, 1])\n",
    "img = render.render_to_image()\n",
    "open3d.io.write_image(\"models/\" + filename + \"_top.png\", img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vis = open3d.visualization.VisualizerWithEditing()\n",
    "vis.create_window()\n",
    "for c in si_meshes:\n",
    "    vis.add_geometry(c)\n",
    "while not vis.is_full_screen():\n",
    "    vis.poll_events()\n",
    "    vis.update_renderer()\n",
    "vis.destroy_window()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vis.destroy_window()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# combine clouds and meshes\n",
    "MasterSiCloud = open3d.geometry.PointCloud()\n",
    "MasterAgCloud = open3d.geometry.PointCloud()\n",
    "MasterSiMesh = open3d.geometry.TriangleMesh()\n",
    "MasterAgMesh = open3d.geometry.TriangleMesh()\n",
    "for c in si_comps:\n",
    "    MasterSiCloud += c\n",
    "for c in ag_comps:\n",
    "    MasterAgCloud += c\n",
    "for c in si_meshes:\n",
    "    MasterSiMesh += c\n",
    "for c in ag_meshes:\n",
    "    MasterAgMesh += c\n",
    "MasterCloud = MasterSiCloud + MasterAgCloud\n",
    "MasterMesh = MasterSiMesh + MasterAgMesh\n",
    "\n",
    "# save clouds and meshes\n",
    "open3d.io.write_point_cloud('models/' + filename + '.pcd', MasterCloud)\n",
    "open3d.io.write_triangle_mesh('models/' + filename + '.ply', MasterMesh)\n",
    "open3d.io.write_point_cloud('models/' + filename + '_Si.pcd', MasterSiCloud)\n",
    "open3d.io.write_triangle_mesh('models/' + filename + '_Si.ply', MasterSiMesh)\n",
    "open3d.io.write_point_cloud('models/' + filename + '_Ag.pcd', MasterAgCloud)\n",
    "open3d.io.write_triangle_mesh('models/' + filename + '_Ag.ply', MasterAgMesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyVista"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Test entire surface '''\n",
    "\n",
    "grid_fill = scipy.ndimage.binary_fill_holes(np.tile(np.where(grid_full > 0, 1, 0), (1,1,1)))\n",
    "grid_surf = hollow(grid_fill)\n",
    "root = remove_unanchored(grid_surf, grid_fill)\n",
    "sroot = sparseFromMatrix(root)\n",
    "point_cloud = pv.PolyData(sroot[:,0:3])\n",
    "point_cloud['species'] = sroot[:,3]\n",
    "#point_cloud['vectors'] = calc_normals(sroot[:,:3])\n",
    "point_cloud.compute_normals()\n",
    "point_cloud.plot(eye_dome_lighting=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#volume = point_cloud.delaunay_3d(progress_bar=True)\n",
    "#shell = volume.extract_geometry()\n",
    "#shell.plot()\n",
    "\n",
    "volume2 = point_cloud.reconstruct_surface(nbr_sz=10, sample_spacing=None, progress_bar=True)\n",
    "shell2 = volume2.extract_geometry()\n",
    "shell2.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Find connected regions '''\n",
    "\n",
    "print(grid_full.shape)\n",
    "grid_si = np.where(grid_full==1, 1, 0)\n",
    "grid_ag = np.where(grid_full==2, 1, 0)\n",
    "regions1, number_of_components1 = scipy.ndimage.label(np.tile(grid_si, (1, 1, 1)), np.ones((3,3,3)))\n",
    "regions2, number_of_components2 = scipy.ndimage.label(np.tile(grid_ag, (1, 1, 1)), np.ones((3,3,3)))\n",
    "regions2 = np.where(regions2>0, regions2 + number_of_components1, 0)\n",
    "regions = regions1 + regions2\n",
    "print(regions1.shape)\n",
    "#print(regions1, number_of_components)\n",
    "\n",
    "comps = np.zeros((grid.shape[0]))\n",
    "for i in range(grid.shape[0]):\n",
    "    p = grid[i,:]\n",
    "    comps[i] = regions[p[2], p[1], p[0]]\n",
    "point_cloud['comps'] = comps\n",
    "point_cloud.plot(eye_dome_lighting=True)\n",
    "\n",
    "\n",
    "#slices = np.zeros(grid_full)\n",
    "#for i in range(grid_full.shape[0]):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct = np.zeros((3,3,3))\n",
    "struct[1,1,:] = 1\n",
    "struct[1,:,1] = 1\n",
    "struct[:,1,1] = 1\n",
    "print(struct)\n",
    "\n",
    "grid_si_fill = scipy.ndimage.binary_fill_holes(np.tile(grid_si, (1,1,1)), structure=struct)\n",
    "grid_ag_fill = scipy.ndimage.binary_fill_holes(np.tile(grid_ag, (1,1,1)), structure=struct)\n",
    "grid_fill = scipy.ndimage.binary_fill_holes(np.tile(np.where(grid_full > 0, 1, 0), (1,1,1)), structure=struct)\n",
    "\n",
    "#grid_si_fill = scipy.ndimage.binary_fill_holes(grid_si_fill, structure=struct)\n",
    "#grid_ag_fill = scipy.ndimage.binary_fill_holes(grid_ag_fill, structure=struct)\n",
    "#grid_fill = scipy.ndimage.binary_fill_holes(grid_fill, structure=struct)\n",
    "\n",
    "#grid_fill = grid_fill[:,grid_fill.shape[1]//3:(grid_fill.shape[1]*2)//3, grid_fill.shape[2]//3:(grid_fill.shape[2]*2)//3]\n",
    "print(grid_fill.shape)\n",
    "\n",
    "grid_interior = copy.deepcopy(grid_fill)\n",
    "for k in range(grid_full.shape[0]):\n",
    "    for j in range(grid_full.shape[1]):\n",
    "        for i in range(grid_full.shape[2]):\n",
    "            if grid_full[k,j,i] != 1:\n",
    "                continue\n",
    "            if grid_full[k,(j-1)%grid_full.shape[1],i] == 0 or grid_full[k,j,(i-1)%grid_full.shape[2]] == 0 or grid_full[k,(j+1)%grid_full.shape[1],i] == 0 or grid_full[k,j,(i+1)%grid_full.shape[2]] == 0:\n",
    "                continue\n",
    "            if k > 0:\n",
    "                if grid_full[k-1,j,i] == 0:\n",
    "                    continue\n",
    "            if k < grid_full.shape[0]-1:\n",
    "                if grid_full[k+1,j,i] == 0:\n",
    "                    continue\n",
    "            grid_interior[k,j,i] = 0\n",
    "hollow = sparseFromMatrix(grid_interior)\n",
    "hollow_cloud = pv.PolyData(hollow[:,0:3])\n",
    "#hollow_cloud['species'] = grid[:,3]\n",
    "hollow_cloud.plot(eye_dome_lighting=True)\n",
    "print(point_cloud.number_of_points, hollow_cloud.number_of_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Remove disconnected components '''\n",
    "\n",
    "regions, number_of_components = scipy.ndimage.label(grid_interior, structure=np.ones((3,3,3)))\n",
    "print(regions.shape)\n",
    "inv_regions = np.where(regions == 1, 1, 0)\n",
    "print(grid_interior.shape, inv_regions.shape)\n",
    "root = grid_interior * inv_regions\n",
    "connected = sparseFromMatrix(root)\n",
    "connected_cloud = pv.PolyData(connected[:,0:3])\n",
    "#connected_cloud['species'] = grid[:,3]\n",
    "connected_cloud.plot(eye_dome_lighting=True)\n",
    "print(point_cloud.number_of_points, connected_cloud.number_of_points)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Ray tracing to find points that can be removed '''\n",
    "root_hollow = copy.deepcopy(root)\n",
    "for k in range(root.shape[0]):\n",
    "    for j in range(root.shape[1]):\n",
    "        for i in range(root.shape[2]):\n",
    "            if(root[k,j,i]) == 0:\n",
    "                continue\n",
    "            # x\n",
    "            if (root[k,j,(i-2)%root.shape[2]] == 0 and root[k,j,(i-1)%root.shape[2]] == 0) or (root[k,j,(i+2)%root.shape[2]] == 0 and root[k,j,(i+1)%root.shape[2]] == 0):\n",
    "                continue\n",
    "            # y\n",
    "            if (root[k,(j-2)%root.shape[1],i] == 0 and root[k,(j-1)%root.shape[1],i] == 0) or (root[k,(j+2)%root.shape[1],1] == 0 and root[k,(j+1)%root.shape[1],i] == 0):\n",
    "                continue\n",
    "            # vertical\n",
    "            if k > 1 and root[k-1,j,i] == 0:\n",
    "                if k > 2 and root[k-2,j,i] == 0:\n",
    "                    continue\n",
    "                continue\n",
    "            # xy, same sign\n",
    "            if (root[k,(j-2)%root.shape[1],(i-2)%root.shape[2]] == 0 and root[k,(j-1)%root.shape[1],(i-1)%root.shape[2]] == 0) or (root[k,(j+2)%root.shape[1],(i+2)%root.shape[2]] == 0 and root[k,(j+1)%root.shape[1],(i+1)%root.shape[2]] == 0):\n",
    "                continue\n",
    "            # xy, diff sign\n",
    "            if (root[k,(j-2)%root.shape[1],(i+2)%root.shape[2]] == 0 and root[k,(j-1)%root.shape[1],(i+1)%root.shape[2]] == 0) or (root[k,(j+2)%root.shape[1],(i-2)%root.shape[2]] == 0 and root[k,(j+1)%root.shape[1],(i-1)%root.shape[2]] == 0):\n",
    "                continue\n",
    "            # diagonals with attitude (do later because annoying)\n",
    "            #if k > 1\n",
    "\n",
    "            # this point cannot be seen so just remove it\n",
    "            root_hollow[k,j,i] = 0\n",
    "    \n",
    "los = sparseFromMatrix(root_hollow)\n",
    "los_cloud = pv.PolyData(los[:,0:3])\n",
    "#los_cloud['species'] = grid[:,3]\n",
    "los_cloud.plot(eye_dome_lighting=True)\n",
    "print(point_cloud.number_of_points, los_cloud.number_of_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Try to create a mesh '''\n",
    "#volume = los_cloud.delaunay_3d(progress_bar=True)\n",
    "#shell = volume.extract_geometry()\n",
    "#shell.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' Maybe a surface from the points? '''\n",
    "#volume2 = los_cloud.reconstruct_surface(nbr_sz=10, sample_spacing=None, progress_bar=True)\n",
    "#shell2 = volume2.extract_geometry()\n",
    "#shell2.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyvista",
   "language": "python",
   "name": "pyvista"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}