{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Marching Cubes {#marching_cubes_example}\n\nGenerate a surface from a scalar field using the flying edges and\nmarching cubes filters as provided by the `contour\n<pyvista.DataSetFilters.contour>`{.interpreted-text role=\"func\"} filter.\n\nSpecial thanks to GitHub user [stla](https://gist.github.com/stla) for\nproviding examples.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nimport pyvista as pv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Spider Cage\n\nUse the marching cubes algorithm to extract the isosurface generated\nfrom the spider cage function.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = 0.9\n\n\ndef spider_cage(x, y, z):\n    x2 = x * x\n    y2 = y * y\n    x2_y2 = x2 + y2\n    return (np.sqrt((x2 - y2) ** 2 / x2_y2 + 3 * (z * np.sin(a)) ** 2) - 3) ** 2 + 6 * (\n        np.sqrt((x * y) ** 2 / x2_y2 + (z * np.cos(a)) ** 2) - 1.5\n    ) ** 2\n\n\n# create a uniform grid to sample the function with\nn = 100\nx_min, y_min, z_min = -5, -5, -3\ngrid = pv.ImageData(\n    dimensions=(n, n, n),\n    spacing=(abs(x_min) / n * 2, abs(y_min) / n * 2, abs(z_min) / n * 2),\n    origin=(x_min, y_min, z_min),\n)\nx, y, z = grid.points.T\n\n# sample and plot\nvalues = spider_cage(x, y, z)\nmesh = grid.contour([1], values, method='marching_cubes')\ndist = np.linalg.norm(mesh.points, axis=1)\nmesh.plot(scalars=dist, smooth_shading=True, cmap=\"plasma\", show_scalar_bar=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Barth Sextic\n\nUse the flying edges algorithm to extract the isosurface generated from\nthe Barth sextic function.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "phi = (1 + np.sqrt(5)) / 2\nphi2 = phi * phi\n\n\ndef barth_sextic(x, y, z):\n    x2 = x * x\n    y2 = y * y\n    z2 = z * z\n    arr = (\n        3 * (phi2 * x2 - y2) * (phi2 * y2 - z2) * (phi2 * z2 - x2)\n        - (1 + 2 * phi) * (x2 + y2 + z2 - 1) ** 2\n    )\n    nan_mask = x2 + y2 + z2 > 3.1\n    arr[nan_mask] = np.nan\n    return arr\n\n\n# create a uniform grid to sample the function with\nn = 100\nk = 2.0\nx_min, y_min, z_min = -k, -k, -k\ngrid = pv.ImageData(\n    dimensions=(n, n, n),\n    spacing=(abs(x_min) / n * 2, abs(y_min) / n * 2, abs(z_min) / n * 2),\n    origin=(x_min, y_min, z_min),\n)\nx, y, z = grid.points.T\n\n# sample and plot\nvalues = barth_sextic(x, y, z)\nmesh = grid.contour([0], values, method='flying_edges')\ndist = np.linalg.norm(mesh.points, axis=1)\nmesh.plot(scalars=dist, smooth_shading=True, cmap=\"plasma\", show_scalar_bar=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Animate Barth Sextic\n\nShow 20 frames of various isocurves extracted from the Barth sextic\nfunction.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def angle_to_range(angle):\n    return -2 * np.sin(angle)\n\n\npl = pv.Plotter(window_size=[800, 800], off_screen=True)\n\npl.open_gif('barth_sextic.gif')\n\nfor angle in np.linspace(0, np.pi, 20, endpoint=False):\n    # clear the plotter before adding each frame's mesh\n    pl.clear()\n    pl.enable_lightkit()\n    mesh = grid.contour([angle_to_range(angle)], values, method='flying_edges')\n    dist = np.linalg.norm(mesh.points, axis=1)\n    pl.add_mesh(\n        mesh,\n        scalars=dist,\n        smooth_shading=True,\n        rng=[0.5, 1.5],\n        cmap=\"plasma\",\n        show_scalar_bar=False,\n    )\n    pl.write_frame()\n\npl.close()"
      ]
    }
  ],
  "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.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}