Merge pull request #90 from autotwin/scale

Scale

.github/workflows/python.yml

@@ -35,7 +35,7 @@ jobs:
     - name: pycodestyle
       run: pycodestyle --verbose .
     - name: pytest
-      run: pytest --verbose .
+      run: pytest --verbose
   wheels:
     strategy:
       fail-fast: false

book/SUMMARY.md

@@ -5,16 +5,20 @@
 # User Guide
 
 - [Installation](installation.md)
-- [Mathematics](mathematics.md)
 - [Examples](examples/README.md)
   - [Unit Tests](examples/unit_tests/README.md)
   - [Spheres](examples/spheres/README.md)
+  - [Spheres - continued](examples/spheres_cont/README.md)
 
 # Reference Guide
 
-- [Command Line Interface](cli.md)
-- [Python Interface](python.md)
-- [Rust Interface](rust.md)
+- [Interfaces](interfaces/README.md)
+  - [Command Line Interface](interfaces/cli.md)
+  - [Python Interface](interfaces/python.md)
+  - [Rust Interface](interfaces/rust.md)
+- [Theory](theory/README.md)
+  - [Mathematics](theory/mathematics.md)
+  - [Smoothing](theory/smoothing.md) 
 
 # Appendices
 

book/examples/spheres/README.md

@@ -15,6 +15,8 @@ resolution (`radius=1`, `radius=3`, and `radius=5`), as shown below:
 
 ![spheres.png](spheres.png)
 
+Figure: Sphere segmentations at selected resolutions, shown in the voxel domain.
+
 For the `radius=1` case, the underyling data structure appears as:
 
 ```python
@@ -95,8 +97,7 @@ array([[[0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0]]], dtype=uint8)
 ```
 
-Because of the large size of `sphere_5`, its data structure is not shown
-here.
+Because of its large size, the data structure for `sphere_5` is not shown here.
 
 These segmentations are saved to
 
@@ -134,5 +135,4 @@ The `spheres_radius_3.inp` file:
 <!-- cmdrun cat spheres_radius_3.inp -->
 ```
 
-Because of the large size of `sphere_5`, its mesh structure is not shown
-here.
+Because of its large size, the mesh structure for `sphere_5` is not shown here.

book/examples/spheres/spheres.py

@@ -1,4 +1,9 @@
 """This module creates a voxelized sphere and exports it as a .npy file.
+
+Example
+-------
+source ~/autotwin/automesh/.venv/bin/activate
+python spheres.py
 """
 
 from pathlib import Path
@@ -95,33 +100,33 @@ def sphere(radius: int, dtype=np.uint8) -> np.ndarray:
 # User input end
 
 
+N_SUBPLOTS = len(spheres)
 IDX = 1
-for title, struc in spheres.items():
-    # for index (key, value) in enumerate(spheres.items()):
-    ax = fig.add_subplot(1, 3, IDX, projection=Axes3D.name)
+for index, (key, value) in enumerate(spheres.items()):
+    ax = fig.add_subplot(1, N_SUBPLOTS, index+1, projection=Axes3D.name)
     ax.voxels(
-        struc,
-        facecolors=colors[IDX-1],
-        edgecolor=colors[IDX-1],
+        value,
+        facecolors=colors[index],
+        edgecolor=colors[index],
         alpha=voxel_alpha,
         lightsource=lightsource)
-    ax.set_title(title)
+    ax.set_title(key)
     IDX += 1
 
     # Set labels for the axes
-    ax.set_xlabel("x")
-    ax.set_ylabel("y")
-    ax.set_zlabel("z")
+    ax.set_xlabel("x (voxels)")
+    ax.set_ylabel("y (voxels)")
+    ax.set_zlabel("z (voxels)")
 
     # Set the camera view
     ax.set_aspect("equal")
     ax.view_init(elev=el, azim=az, roll=roll)
 
     if serialize:
-        cc = aa.with_stem("spheres_" + title)
+        cc = aa.with_stem("spheres_" + key)
         dd = cc.with_suffix(".npy")
         # Save the data in .npy format
-        np.save(dd, struc)
+        np.save(dd, value)
         print(f"Saved: {dd}")
 
 fig.tight_layout()

book/examples/spheres_cont/README.md

@@ -0,0 +1,66 @@
+# Spheres - Continued
+
+Use the fundamentals learned in the [previous example](../spheres/README.md) to create a more sophisticated example:  Concentric, high-resolution spheres consisting of three materials.
+
+## Problem Statement
+
+### Given
+
+Given three concentric spheres of radius 10, 11, and 12 cm, as shown in the figure below,
+
+![spheres_cont_dim](spheres_cont_dim.png)
+
+Figure: Schematic cross-section of three concentric spheres of radius 10, 11, and 12 cm.  Grid spacing is 1 cm.
+
+### Find
+
+Use the following segmentation resolutions,
+
+resolution (vox/cm) | element side length (cm) | `nelx` | # voxels
+---: | :---: | ---: | ---:
+1 | 1.0 | 24 | 13,824
+2 | 0.5 | 48 | 110,592
+4 | 025 | 96 | 884,736
+10 | 0.1 | 240 | 13,824,000
+
+with a cubic domain (`nelx = nely = nelz`),
+to create finite element meshes.
+
+## Solution
+
+### Python Segmentation
+
+Use [spheres_cont.py](spheres_cont.py) to create segmentations,
+
+```python
+<!-- cmdrun cat spheres_cont.py -->
+```
+
+![spheres_cont](spheres_cont.png)
+
+Figure: Sphere segmentations at selected resolutions, shown in the voxel domain.
+
+### Autotwin
+
+Use `Autotwin` to convert the segmentations into finite element meshes.
+
+```sh
+automesh -i spheres_resolution_1.npy -o spheres_resolution_1.inp -x 24 -y 24 -z 24 -xtranslate -12 -ytranslate -12 -ztranslate -12
+# didn't work
+
+automesh -i spheres_resolution_1.npy -o spheres_resolution_1.inp -x 24 -y 24 -z 24 --xtranslate 24 --ytranslate 24 --ztranslate 24
+# works
+```
+
+But translating in a negative direction seems to be error prone:
+
+```sh
+automesh -i spheres_resolution_1.npy -o spheres_resolution_1.inp -x 24 -y 24 -z 24 --xtranslate -24 --ytranslate 24 --ztranslate 24
+# didn't work
+
+error: unexpected argument '-2' found
+```
+
+### Meshes
+
+Figure (to come): Sphere segmentations at selected resolutions, shown in the finite element domain with cutting plane to expose interior structure.

book/examples/spheres_cont/spheres_cont.py

@@ -0,0 +1,165 @@
+"""This module builds on the `spheres.py` module to create high resolution,
+three-material, concentric spheres and export the voxelization as a .npy
+file.
+
+Example
+-------
+source ~/autotwin/automesh/.venv/bin/activate
+python spheres_cont.py
+"""
+
+from pathlib import Path
+from typing import Final
+
+from matplotlib.colors import LightSource
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import numpy as np
+
+
+def sphere(resolution: int, dtype=np.uint8) -> np.ndarray:
+    """Generate a 3D voxelized representation of three concentric spheres
+        of 10, 11, and 12 cm, at a given resolution.
+
+        Parameters
+        ----------
+        resolution : int
+            The resolution as voxels per centimeter.  Minimum value is 1.
+
+        dtype: data-type, optional
+            The data type of the output array.  Default is np.uint8.
+
+        Returns
+        -------
+        np.ndarray
+            A 3D numpy array representing the voxelized spheres.  Voxels within
+            the inner sphere are set to 1, the intermediate shell are set to 2,
+            and the outer shell are set to 3.  Voxels outside the spheres are
+            set to 0.
+
+        Raises
+        ------
+        ValueError
+            If the resolution is less than 1.
+
+        Example
+        -------
+        >>> sphere(resolution=2) returns
+            array(
+                [
+    #                [[0, 0, 0], [0, 1, 0], [0, 0, 0]],
+    #                [[0, 1, 0], [1, 1, 1], [0, 1, 0]],
+    #                [[0, 0, 0], [0, 1, 0], [0, 0, 0]]
+                ],
+                dtype=uint8
+            )
+    """
+    if resolution < 1:
+        raise ValueError("Resolution must be >= 1")
+
+    r10 = 10  # cm
+    r11 = 11  # cm
+    r12 = 12  # cm
+
+    # We change the algorithm a bit here so we can exactly match the radius:
+    # number of voxels per side length (nvps)
+    # nvps = 2 * r12 * resolution + 1
+    nvps = 2 * r12 * resolution
+    vox_z, vox_y, vox_x = np.mgrid[
+        -r12: r12: nvps * 1j,
+        -r12: r12: nvps * 1j,
+        -r12: r12: nvps * 1j,
+    ]
+    domain = vox_x**2 + vox_y**2 + vox_z**2
+    mask_10_in = np.array(domain <= r10 * r10, dtype=dtype)
+    mask_11_in = np.array(domain <= r11 * r11, dtype=dtype)
+    mask_12_in = np.array(domain <= r12 * r12, dtype=dtype)
+
+    mask_10_11 = mask_11_in - mask_10_in
+    mask_11_12 = mask_12_in - mask_11_in
+
+    shell_10_11 = 2 * mask_10_11
+    shell_11_12 = 3 * mask_11_12
+
+    result = mask_10_in + shell_10_11 + shell_11_12
+    # breakpoint()
+    return result
+
+
+rr = (1, 2, 4, 10)  # resolutions (voxels per cm)
+lims = tuple(map(lambda x: [0, 24*x], rr))  # limits
+tt = tuple(map(lambda x: [0, 12*x, 24*x], rr))  # ticks
+
+# User input begin
+spheres = {
+    "resolution_1": sphere(resolution=rr[0]),
+    "resolution_2": sphere(resolution=rr[1]),
+}
+
+aa = Path(__file__)
+bb = aa.with_suffix(".png")
+
+# Visualize the elements.
+# width, height = 8, 4
+width, height = 6, 3
+fig = plt.figure(figsize=(width, height))
+# fig = plt.figure(figsize=(8, 8))
+
+el, az, roll = 63, -110, 0
+cmap = plt.get_cmap(name="tab10")
+num_colors = len(spheres)
+voxel_alpha: Final[float] = 0.9
+
+colors = cmap(np.linspace(0, 1, num_colors))
+lightsource = LightSource(azdeg=325, altdeg=45)  # azimuth, elevation
+# lightsource = LightSource(azdeg=325, altdeg=90)  # azimuth, elevation
+dpi: Final[int] = 300  # resolution, dots per inch
+visualize: Final[bool] = True  # turn to True to show the figure on screen
+serialize: Final[bool] = True  # turn to True to save .png and .npy files
+# User input end
+
+N_SUBPLOTS = len(spheres)
+for index, (key, value) in enumerate(spheres.items()):
+    print(f"index: {index}")
+    print(f"key: {key}")
+    print(f"value: {value}")
+    ax = fig.add_subplot(1, N_SUBPLOTS, index+1, projection=Axes3D.name)
+    ax.voxels(
+        value,
+        facecolors=colors[index],
+        edgecolor=colors[index],
+        alpha=voxel_alpha,
+        lightsource=lightsource)
+    ax.set_title(key)
+
+    # Set labels for the axes
+    ax.set_xlabel("x (voxels)")
+    ax.set_ylabel("y (voxels)")
+    ax.set_zlabel("z (voxels)")
+
+    ax.set_xticks(ticks=tt[index])
+    ax.set_yticks(ticks=tt[index])
+    ax.set_zticks(ticks=tt[index])
+
+    ax.set_xlim(lims[index])
+    ax.set_ylim(lims[index])
+    ax.set_zlim(lims[index])
+
+    # Set the camera view
+    ax.set_aspect("equal")
+    ax.view_init(elev=el, azim=az, roll=roll)
+
+    if serialize:
+        cc = aa.with_stem("spheres_" + key)
+        dd = cc.with_suffix(".npy")
+        # Save the data in .npy format
+        np.save(dd, value)
+        print(f"Saved: {dd}")
+
+# fig.tight_layout()  # don't use as it clips the x-axis label
+if visualize:
+    plt.show()
+
+if serialize:
+    fig.savefig(bb, dpi=dpi)
+    print(f"Saved: {bb}")

book/examples/unit_tests/README.md

@@ -1,6 +1,11 @@
 # Tests
 
-Following is documentation for unit tests used to validate code implementation.
+The following illustates a subset of the unit tests used to validate
+the code implementation.
+For a complete listing of the unit sets, see
+[voxels.rs](https://github.com/autotwin/automesh/blob/main/tests/voxel.rs)
+and [voxel.py](https://github.com/autotwin/automesh/blob/main/tests/voxel.py).
+
 The Python script used to generate the figures is included [below](#source).
 
 **Remark:** We use the convention `np` when importing `numpy` as follows:

book/interfaces/README.md

@@ -0,0 +1,7 @@
+# Interfaces
+
+There are three main ways to use `autotwin`:
+
+* The [command line interface](cli.md)
+* The [Python interface](python.md)
+* The [Rust interface](rust.md)

book/theory/README.md

@@ -0,0 +1,4 @@
+# Theory
+
+This section contains a description of some of the theory used
+for implementation.