Merge pull request #83 from autotwin/spheres

Spheres

book/examples/spheres/README.md

@@ -1,8 +1,10 @@
 # Spheres
 
+We segment a sphere into coarse voxel meshes.
+
 ## Segmentation
 
-Using [spheres.py](spheres.py), 
+Using [spheres.py](spheres.py),
 
 ```python
 <!-- cmdrun cat spheres.py -->
@@ -16,7 +18,7 @@ resolution (`radius=1`, `radius=3`, and `radius=5`), as shown below:
 For the `radius=1` case, the underyling data structure appears as:
 
 ```python
-spheres["ball_1"]
+spheres["sphere_1"]
 
 array([[[0, 0, 0],
         [0, 1, 0],
@@ -34,7 +36,7 @@ array([[[0, 0, 0],
 For the `radius=3` case, the underyling data structure appears as:
 
 ```python
-spheres["ball_3"]
+spheres["sphere_3"]
 
 array([[[0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0],
@@ -93,25 +95,45 @@ array([[[0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0]]], dtype=uint8)
 ```
 
-Because of the large size of `ball_5`, its data structure is not shown
+Because of the large size of `sphere_5`, its data structure is not shown
 here.
 
 These data structures are saved to
 
-* [spheres_ball_1_.npy](spheres_ball_1_.npy)
-* [spheres_ball_3_.npy](spheres_ball_3_.npy)
-* [spheres_ball_5_.npy](spheres_ball_5_.npy)
+* [spheres_sphere_1.npy](spheres_sphere_1.npy)
+* [spheres_sphere_3.npy](spheres_sphere_3.npy)
+* [spheres_sphere_5.npy](spheres_sphere_5.npy)
+
+## Autotwin
+
+```sh
+cargo run -- -i book/examples/spheres/spheres_sphere_1.npy -o book/examples/spheres/spheres_sphere_1.inp -x 3 -y 3 -z 3
+```
+
+```sh
+cargo run -- -i book/examples/spheres/spheres_sphere_3.npy -o book/examples/spheres/spheres_sphere_3.inp -x 7 -y 7 -z 7
+```
+
+```sh
+cargo run -- -i book/examples/spheres/spheres_sphere_5.npy -o book/examples/spheres/spheres_sphere_5_.inp -x 11 -y 11 -z 11
+```
 
 ## Mesh
 
+The `spheres_sphere_1.inp` file:
+
 ```sh
-cargo run -- -i book/examples/spheres/spheres_ball_1_.npy -o book/examples/spheres/spheres_ball_1_.inp -x 3 -y 3 -z 3
+<!-- cmdrun cat spheres_sphere_1.inp -->
 ```
 
+The `spheres_sphere_3.inp` file:
+
 ```sh
-cargo run -- -i book/examples/spheres/spheres_ball_3_.npy -o book/examples/spheres/spheres_ball_3_.inp -x 7 -y 7 -z 7
+<!-- cmdrun cat spheres_sphere_3.inp -->
 ```
 
+The `spheres_sphere_5.inp` file:
+
 ```sh
-cargo run -- -i book/examples/spheres/spheres_ball_5_.npy -o book/examples/spheres/spheres_ball_5_.inp -x 11 -y 11 -z 11
+<!-- cmdrun cat spheres_sphere_5.inp -->
 ```

book/examples/spheres/spheres.py

@@ -1,9 +1,5 @@
 """This module demonstrates how to create a voxelized sphere and export
 it as a .npy file.
-
-Items that need to be installed into the virtual environment:
-    pip install scikit-image
-
 """
 
 import matplotlib.pyplot as plt
@@ -13,15 +9,67 @@
 from pathlib import Path
 from typing import Final
 
-from skimage.morphology import ball
+
+def sphere(radius: int, dtype=np.uint8) -> np.ndarray:
+    """Generate a 3D voxelized representation of a sphere.
+
+    Parameters
+    ----------
+    radius: int
+        The radius of the sphere.  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 of shape (2*radius+1, 2*radius+1, 2*radius+1)
+        representing the voxelized sphere.  Voxels within the sphere are
+        set to 1, and those outside are set to 0.
+
+    Raises
+    ------
+    ValueError
+        If the radius is less than 1.
+
+    Example
+    -------
+    >>> sphere(radius=1) 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
+        )
+
+    Reference
+    ---------
+    Adapted from:
+    https://github.com/scikit-image/scikit-image/blob/v0.24.0/skimage/morphology/footprints.py#L763-L833
+    """
+    if radius < 1:
+        raise ValueError("Radius must be >= 1")
+
+    n_voxels_per_side = 2 * radius + 1
+    vox_z, vox_y, vox_x = np.mgrid[
+        -radius:radius:n_voxels_per_side * 1j,
+        -radius:radius:n_voxels_per_side * 1j,
+        -radius:radius:n_voxels_per_side * 1j,
+    ]
+    voxel_radius_squared = vox_x**2 + vox_y**2 + vox_z**2
+    result = np.array(voxel_radius_squared <= radius * radius, dtype=dtype)
+    return result
 
 
 # User input begin
 
 spheres = {
-    "ball_1": ball(radius=1),
-    "ball_3": ball(radius=3),
-    "ball_5": ball(radius=5),
+    "sphere_1": sphere(radius=1),
+    "sphere_3": sphere(radius=3),
+    "sphere_5": sphere(radius=5),
 }
 
 aa = Path(__file__)
@@ -42,12 +90,14 @@
 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] = False  # turn to True to show the figure on screen
 serialize: Final[bool] = False  # turn to True to save .png and .npy files
 # User input end
 
 
 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)
     ax.voxels(
         struc,
@@ -68,14 +118,15 @@
     ax.view_init(elev=el, azim=az, roll=roll)
 
     if serialize:
-        cc = aa.with_stem("spheres_" + title + "_")
+        cc = aa.with_stem("spheres_" + title)
         dd = cc.with_suffix(".npy")
         # Save the data in .npy format
         np.save(dd, struc)
         print(f"Saved: {dd}")
 
 fig.tight_layout()
-plt.show()
+if visualize:
+    plt.show()
 
 if serialize:
     fig.savefig(bb, dpi=dpi)

book/examples/spheres/spheres_ball_1_.inp → book/examples/spheres/spheres_sphere_1.inp

@@ -1,7 +1,7 @@
 *HEADING
 autotwin.automesh
-version 0.1.3
-autogenerated on 2024-09-11 23:05:42.323084 UTC
+version 0.1.6
+autogenerated on 2024-09-19 23:00:58.449967 UTC
 **
 *NODE, NSET=ALLNODES
      1,      1.000000e0,      1.000000e0,      0.000000e0

book/examples/spheres/spheres_ball_3_.inp → book/examples/spheres/spheres_sphere_3.inp

@@ -1,7 +1,7 @@
 *HEADING
 autotwin.automesh
-version 0.1.3
-autogenerated on 2024-09-11 23:05:49.991723 UTC
+version 0.1.6
+autogenerated on 2024-09-19 23:01:48.003423 UTC
 **
 *NODE, NSET=ALLNODES
       1,      3.000000e0,      3.000000e0,      0.000000e0

book/examples/spheres/spheres_ball_5_.inp → book/examples/spheres/spheres_sphere_5.inp

@@ -1,7 +1,7 @@
 *HEADING
 autotwin.automesh
-version 0.1.3
-autogenerated on 2024-09-11 23:05:56.317604 UTC
+version 0.1.6
+autogenerated on 2024-09-19 23:02:10.991462 UTC
 **
 *NODE, NSET=ALLNODES
       1,      5.000000e0,      5.000000e0,      0.000000e0

book/examples/spheres/test_spheres.py

@@ -0,0 +1,152 @@
+"""This module performs point testing of the sphere module.
+
+Example
+-------
+
+> python -m pytest book/examples/spheres/test_spheres.py
+
+"""
+
+import numpy as np
+import pytest
+
+import spheres as sph
+
+
+def test_sphere():
+    """Unit tests for the sphere function."""
+
+    # Assure that radius >=1 assert is raised
+    with pytest.raises(ValueError, match="Radius must be >= 1"):
+        sph.sphere(radius=0)
+
+    #  Assure radius=1 is correct
+    gold_r1 = np.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=np.uint8,
+    )
+
+    result_r1 = sph.sphere(radius=1)
+    assert np.all(gold_r1 == result_r1)
+
+    #  Assure radius=2 is correct
+    gold_r2 = np.array(
+        [
+            [
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 1, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0],
+                [0, 1, 1, 1, 0],
+                [0, 1, 1, 1, 0],
+                [0, 1, 1, 1, 0],
+                [0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 1, 0, 0],
+                [0, 1, 1, 1, 0],
+                [1, 1, 1, 1, 1],
+                [0, 1, 1, 1, 0],
+                [0, 0, 1, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0],
+                [0, 1, 1, 1, 0],
+                [0, 1, 1, 1, 0],
+                [0, 1, 1, 1, 0],
+                [0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 1, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+            ],
+        ],
+        dtype=np.uint8,
+    )
+
+    result_r2 = sph.sphere(radius=2)
+    assert np.all(gold_r2 == result_r2)
+
+    #  Assure radius=3 is correct
+    gold_r3 = np.array(
+        [
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 1, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 1, 1, 1, 0, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 0, 1, 1, 1, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 1, 0, 0, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [1, 1, 1, 1, 1, 1, 1],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 0, 0, 1, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 1, 1, 1, 0, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 1, 1, 1, 1, 1, 0],
+                [0, 0, 1, 1, 1, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+            [
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 1, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0],
+            ],
+        ],
+        dtype=np.uint8,
+    )
+
+    result_r3 = sph.sphere(radius=3)
+    assert np.all(gold_r3 == result_r3)

pyproject.toml

@@ -28,7 +28,6 @@ dev = [
   'pycodestyle',
   'pytest',
   'pytest-cov',
-  'scikit-image',
 ]
 
 [project.urls]