--- title: "Draw an isosurface with Python Plotly with the help of Pyvista" author: "Stéphane Laurent" date: '2022-01-25' output: md_document: variant: markdown preserve_yaml: yes html_document: highlight: kate keep_md: no tags: geometry, plotly, graphics, python, pyvista highlighter: pandoc-solarized --- **Plotly** isosurfaces are relatively new. I have not tried yet in Python. I tried in R and I found them a bit slow. [In this post](https://laustep.github.io/stlahblog/posts/plotly_trisurf.html) I showed how to plot an isosurface in R **Plotly** with the help of the **misc3d** package: by constructing a triangular 3D mesh and using the `mesh3d` type in **Plotly**. Here I will show how to do something similar in Python, with the help of the **PyVista** package. I don't know whether this is faster than the built-in way to draw an isosurface in **Plotly**, but there's at least one advantage: the one I mentioned [here](https://laustep.github.io/stlahblog/posts/MeshClipping.html), that is, the possibility to clip the isosurface to a region (usually a ball). I don't think this is possible with the **Plotly** way (I may be wrong). Let's go. I gonna draw the Togliatti surface again. First, the mesh clipped to a box: ```python from math import sqrt import numpy as np import pyvista as pv import plotly.graph_objects as go import plotly.io as pio # to save the graphics as html def f(x, y, z): return ( 64 * (x - 1) * ( x**4 - 4 * x**3 - 10 * x**2 * y**2 - 4 * x**2 + 16 * x - 20 * x * y**2 + 5 * y**4 + 16 - 20 * y**2 ) - 5 * sqrt(5 - sqrt(5)) * (2 * z - sqrt(5 - sqrt(5))) * (4 * (x**2 + y**2 - z**2) + (1 + 3 * sqrt(5)))**2 ) # generate data grid for computing the values X, Y, Z = np.mgrid[(-5):5:250j, (-5):5:250j, (-4):4:250j] # create a structured grid grid = pv.StructuredGrid(X, Y, Z) # compute and assign the values values = f(X, Y, Z) grid.point_data["values"] = values.ravel(order = "F") # compute the isosurface f(x, y, z) = 0 isosurf = grid.contour(isosurfaces = [0]) mesh0 = isosurf.extract_geometry() ``` Now we clip the Togliatti surface to a ball. The mesh we obtain is named `mesh`. It is triangular: `mesh.is_all_triangles` returns `True`. If it were not, we would need to run `mesh.triangulate`. Finally we extract the triangular faces from this mesh. The easiest way is to reshape the vector `mesh.faces` to a matrix. ```python # surface clipped to the ball of radius 4.8, with the help of `clip_scalar`: mesh0["dist"] = np.linalg.norm(mesh0.points, axis = 1) mesh = mesh0.clip_scalar("dist", value = 4.8) mesh.is_all_triangles # True # mesh.triangulate() points = mesh.points triangles = mesh.faces.reshape(-1, 4) ``` Now it's a child game to plot the isosurface with **Plotly**: ```python fig = go.Figure(data=[ go.Mesh3d( x=points[:, 0], y=points[:, 1], z=points[:, 2], colorscale = [[0, 'gold'], [0.5, 'mediumturquoise'], [1, 'magenta']], # Intensity of each vertex, which will be interpolated and color-coded intensity = np.linspace(0, 1, len(mesh.points)), # i, j and k give the vertices of the triangles i = triangles[:, 1], j = triangles[:, 2], k = triangles[:, 3], showscale = False ) ]) fig.show() ```