Expanding documentation

docs/background.rst

@@ -0,0 +1,53 @@
+.. _background:
+.. index:: Background
+
+Background
+==========
+
+Constructing an isosurface of a scalar field essentially constitutes finding those
+set of point on the grid that share an **isovalue** and connecting those points.
+Many explanations have been written on the topic and I warmly invite the reader
+to study the excellent
+`work of Paul Bourke <https://paulbourke.net/geometry/polygonise/>`_ on this topic.
+
+Algorithm
+---------
+
+:program:`PyTessel` uses the
+`marching cubes <https://en.wikipedia.org/wiki/Marching_cubes>`_ algorithm.
+In this algorithm,
+the scalar field is represented as a set of point on a (typically regular)
+grid. Each octet of grid points thus represents a cube. For each of the 8
+vertices in each cell, it is determined whether the value of the function is
+larger or smaller than the isovalue. If for one vertex the value is larger
+than the isovalue and for another vertex the value smaller than the isovalue,
+then the isosurface has to cut the edge that is shared by the two vertices.
+
+Given eight vertices per cube, there are only :math:`2^{8}=256` possibilities of
+how a scalar field can interact with each cube. An overview of all these
+possibilities is provided in the image below. Note that within this set of 256
+possibilities, there are only 15 unique results as the majority of the results
+are linked by a symmetry operation such as a rotation. All these 256 different
+ways by which the scalar field can interact with each cell are stored in a
+pre-calculated table by which the resulting list of edge intersections and the
+polygons that originate from these intersections can be readily found. The
+exact position where the isosurface intersects each edge is determined using
+trilinear interpolation.
+
+.. image:: _static/img/marching_cubes.png
+
+The algorithm essentially proceeds as follows:
+
+1. The space wherein the function is evaluated is discretized into small
+   cubes.
+2. For each vertex on these cubes the function is evaluated and it is
+   determined whether the value for the function at each of the vertices is larger
+   or smaller than the isovalue. This leads to a number of edge intersections for
+   each cube from which the polygons (triangles) for each cube can be established.
+3. All polygons are gathered to form the threedimensional isosurface.
+
+It should be noted that this algorithm can be executed in a highly efficient
+fashion using trivial parallellization as the result for each cube is
+completely independent from all the other cubes. It turns out that the
+generation of the scalar field, which by itself is also a highly
+parallellizable step, is typically the most time-consuming.

docs/examples.rst

@@ -0,0 +1,132 @@
+.. _examples:
+.. index:: Examples
+
+Examples
+========
+
+Icosahedral metaballs
+---------------------
+
+Here, we provide a simple example how to use :program:`PyTessel`. In this
+example, we will be constructing the isosurface of a scalar field corresponding
+to 12 `metaballs <https://en.wikipedia.org/wiki/Metaballs>`_ arranged such
+that they lie on the vertices of a regular
+`icosahedron <https://en.wikipedia.org/wiki/Icosahedron>`_.
+
+In this script, we will be exploring the quality of the final isosurface as
+function of the number of sampling points originally taken when constructing
+the underlying scalar field. As can be seen in the :code:`main()` function
+of the script, we loop over the array :code:`[10,20,25,50,100,200]` and execute
+the :code:`marching_cubes` algorithm for differently sampled scalar fields and
+construct a separate isosurface for each scalar field.
+
+.. code:: python
+
+    from pytessel import PyTessel
+    import numpy as np
+
+    def main():
+        """
+        Build 6 .ply files of increasing quality
+        """
+        pytessel = PyTessel()
+
+        for nrpoints in [10,20,25,50,100,200]:
+            sz = 3
+
+            x = np.linspace(-sz, sz, nrpoints)
+            y = np.linspace(-sz, sz, nrpoints)
+            z = np.linspace(-sz, sz, nrpoints)
+
+            xx, yy, zz, field = icosahedron_field(x,y,z)
+
+            unitcell = np.diag(np.ones(3) * sz * 2)
+            pytessel = PyTessel()
+            isovalue = 3.75
+            vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue)
+
+            pytessel.write_ply('icosahedron_%03i.ply' % nrpoints, vertices, normals, indices)
+
+    def icosahedron_field(x,y,z):
+        """
+        Produce a scalar field for the icosahedral metaballs
+        """
+        phi = (1 + np.sqrt(5)) / 2
+        vertices = [
+            [0,1,phi],
+            [0,-1,-phi],
+            [0,1,-phi],
+            [0,-1,phi],
+            [1,phi,0],
+            [-1,-phi,0],
+            [1,-phi,0],
+            [-1,phi,0],
+            [phi,0,1],
+            [-phi,0,-1],
+            [phi,0,-1],
+            [-phi,0,1]
+        ]
+
+        xx,yy,zz = np.meshgrid(x,y,z)
+        field = np.zeros_like(xx)
+        for v in vertices:
+            field += f(xx,yy,zz,v[0], v[1],v[2])
+
+        return xx,yy,zz,field
+
+    def f(x,y,z,X0,Y0,Z0):
+        """
+        Single metaball function
+        """
+        return 1 / ((x - X0)**2 + (y - Y0)**2 + (z - Z0)**2)
+
+    if __name__ == '__main__':
+        main()
+
+Each of the isosurfaces are rendered using
+`Blender <https://www.blender.org/>`_. The result is found in the image below.
+From this image, we can see that upon increasing the number of sample points,
+we gradually increase the quality of the isosurface.
+
+.. image:: _static/img/metaballs_3x2.png
+
+Execution times
+***************
+
+To get an impression on the performance of the program, we here compare the
+execution times for scalar field generation and isosurface constructing as function
+of the number of sampling points. As is typically observed, the generation of
+the scalar field takes more time than the actual construction of the isosurface.
+
+.. list-table:: Execution times for scalar field and isosurface generation.
+   :widths: 25 25 25
+   :header-rows: 1
+
+   * - Number of points
+     - Scalar field generation (seconds)
+     - Isosurface construction (seconds)
+   * - 10
+     - 0.000000
+     - 0.001002
+   * - 20
+     - 0.000499
+     - 0.004001
+   * - 25
+     - 0.000500
+     - 0.006998
+   * - 50
+     - 0.026501
+     - 0.032499
+   * - 100
+     - 0.212501
+     - 0.155500
+   * - 200
+     - 1.865499
+     - 0.817000
+
+We can also assess the complexity of the algorithm by sampling the execution time
+as function of the number of data points as shown in the image below. On the
+basis of the slope, we find that scalar field generation scales by
+:math:`\sim N^{3.5}` whereas isosurface generation scales by :math:`\sim N^{2.5}`.
+
+.. image:: _static/img/scaling_relation.jpg

docs/index.rst

@@ -12,7 +12,13 @@ PyTessel: isosurface generation tool
 
 :program:`PyTessel` is a Python package for building isosurfaces from 3D scalar
 fields. It is originally designed to create isosurfaces of (molecular) orbitals
-and electron densities, but is agnotistic with respect to data input.
+and electron densities, but is agnotistic with respect to data input. The
+majority of the instructions of :program:`PyTessel` are coded in C++ and are
+coupled to a number of Python routines using the `Cython <https://cython.org/>`_
+extension. :program:`PyTessel` makes use of shared-memory parallelization
+via `OpenMP <https://www.openmp.org/>`_ ensuring efficient use is made of
+trivial parallelization strategies. For more information about the algorithm,
+please consult the :ref:`background` section.
 
 To construct an isosurface, :program:`PyTessel` requires a scalar field and a
 description of the unit cell wherein the scalar field resides. A very succinct
@@ -61,7 +67,9 @@ requests are ideally submitted via the `github issue tracker
    :caption: Contents:
 
    installation
+   background
    user_interface
+   examples
    community_guidelines
 
 Indices and tables

example/.gitignore

@@ -1 +1,3 @@
 *.ply
+*.jpg
+*.png

example/metaballs_scaling.py

@@ -0,0 +1,104 @@
+from pytessel import PyTessel
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+
+# In this script, we sample the scaling relationship between the number of
+# sample points and the execution time for both isosurface construction as well
+# as scalar field generation
+
+def main():
+    """
+    Build 6 .ply files of increasing quality
+    """
+    pytessel = PyTessel()
+
+    nrpoints = [20,25,50,100,200,500]
+    scalarfield_times = []
+    isosurface_times = []
+
+    for nrpoint in nrpoints:
+        sz = 3
+
+        x = np.linspace(-sz, sz, nrpoint)
+        y = np.linspace(-sz, sz, nrpoint)
+        z = np.linspace(-sz, sz, nrpoint)
+
+
+        start = time.time()
+        xx, yy, zz, field = icosahedron_field(x,y,z)
+        end = time.time()
+        scalarfield_times.append(end - start)
+
+        unitcell = np.diag(np.ones(3) * sz * 2)
+        pytessel = PyTessel()
+        isovalue = 3.75
+
+        start = time.time()
+        vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue)
+        end = time.time()
+        isosurface_times.append(end - start)
+
+
+    fig, ax = plt.subplots(1, 2, dpi=144)
+    ax[0].loglog(nrpoints, scalarfield_times, 'o', color='red', label='Data points')
+    ax[1].loglog(nrpoints, isosurface_times, 'o', color='red', label='Data points')
+
+    ax[0].set_xlabel('Number of data points [-]')
+    ax[0].set_ylabel('Execution time [s]')
+    ax[1].set_xlabel('Number of data points [-]')
+    ax[1].set_ylabel('Execution time [s]')
+
+    # perform linear fit
+    a,b = np.polyfit(np.log(nrpoints), np.log(scalarfield_times), deg=1)
+    x = np.linspace(np.min(nrpoints), np.max(nrpoints), 100)
+    y = np.exp(a * np.log(x) + b)
+    ax[0].plot(x, y, '--', color='black', label='Linear fit: %.2f' % a)
+    ax[0].set_title('Scalar field generation' % a)
+    ax[0].legend()
+
+    a,b = np.polyfit(np.log(nrpoints), np.log(isosurface_times), deg=1)
+    x = np.linspace(np.min(nrpoints), np.max(nrpoints), 100)
+    y = np.exp(a * np.log(x) + b)
+    ax[1].plot(x, y, '--', color='black', label='Linear fit: %.2f' % a)
+    ax[1].set_title('Isosurface construction' % a)
+    ax[1].legend()
+
+    plt.tight_layout()
+    plt.savefig('scaling_relation.jpg')
+
+def icosahedron_field(x,y,z):
+    """
+    Produce a scalar field for the icosahedral metaballs
+    """
+    phi = (1 + np.sqrt(5)) / 2
+    vertices = [
+        [0,1,phi],
+        [0,-1,-phi],
+        [0,1,-phi],
+        [0,-1,phi],
+        [1,phi,0],
+        [-1,-phi,0],
+        [1,-phi,0],
+        [-1,phi,0],
+        [phi,0,1],
+        [-phi,0,-1],
+        [phi,0,-1],
+        [-phi,0,1]
+    ]
+
+    xx,yy,zz = np.meshgrid(x,y,z)
+    field = np.zeros_like(xx)
+    for v in vertices:
+        field += f(xx,yy,zz,v[0], v[1],v[2])
+
+    return xx,yy,zz,field
+
+def f(x,y,z,X0,Y0,Z0):
+    """
+    Single metaball function
+    """
+    return 1 / ((x - X0)**2 + (y - Y0)**2 + (z - Z0)**2)
+
+if __name__ == '__main__':
+    main()