3D Binary Image Analysis

This tutorial covers 3D binary image analysis. The basic procedures are the same as with 2D binary image analysis, so please complete that tutorial before proceeding with this one.

We will analyze a 50x50x50 pixel 3D image.

There are several ways to represent 3D images, but in this tutorial, we'll use cross-sectional images for each layer. This means we'll be working with 50 individual 50x50 pixel images.

These files are stored in the data/ directory, and their filenames reflect the order of their respective cross-sections.

ls data/

0000.png  0007.png  0014.png  0021.png  0028.png  0035.png  0042.png  0049.png
0001.png  0008.png  0015.png  0022.png  0029.png  0036.png  0043.png
0002.png  0009.png  0016.png  0023.png  0030.png  0037.png  0044.png
0003.png  0010.png  0017.png  0024.png  0031.png  0038.png  0045.png
0004.png  0011.png  0018.png  0025.png  0032.png  0039.png  0046.png
0005.png  0012.png  0019.png  0026.png  0033.png  0040.png  0047.png
0006.png  0013.png  0020.png  0027.png  0034.png  0041.png  0048.png

Loading Data Files

Load these images and create a 50x50x50 pixel binary array from the data.

# Import libraries
import imageio
import numpy as np
import matplotlib.pyplot as plt

Before constructing a 3D array, we will load 0000.png and check what kind of image it is.

image_0000 = imageio.v3.imread("data/0000.png")

image_0000.shape, image_0000.dtype, np.min(image_0000), np.max(image_0000)

((50, 50), dtype('uint8'), np.uint8(0), np.uint8(255))

This is a 50x50 array whose elements are between 0 and 255. That is, the image is gray-scaled.

Check the contents of the array.

image_0000

array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ...,   0,   0,   0],
       [255, 255, 255, ...,   0,   0,   0],
       [255, 255, 255, ...,   0,   0,   0]], shape=(50, 50), dtype=uint8)

plt.imshow(image_0000, "gray")

<matplotlib.image.AxesImage at 0x74bf4bf7ce00>

Load the 50 image files in the order of filenames and stack them.

pict = np.stack([
    imageio.v3.imread("data/{:04d}.png".format(n)) > 128
    for n in range(50) 
], axis=0)

pict.shape, pict.dtype

((50, 50, 50), dtype('bool'))

The 50×50×50 array has been completed. Let's visualize it in 3D.

We'll use Plotly for the visualization. We will use plotly.graph_objects and homcloud.plotly_3d modules. See API manual in more details.

You can change the opacity by enabling the commented-out line.

import homcloud.plotly_3d as p3d
import plotly.graph_objects as go
import plotly.io
plotly.io.renderers.default = 'notebook'  # plotly Setting. Choose your favorite renderer.

fig = go.Figure(data=[p3d.Bitmap3d(pict)], layout=dict(scene=p3d.SimpleScene()))
# fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.show()

Computing the Persistence Diagrams

We will compute the persistence diagrams.
We'll use hc.PDList.from_bitmap_levelset and hc.distance_transform.
We can compute the diagrams in the same way as in the 2D case.

import homcloud.interface as hc

hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="binary-3d.pdgm")

PDList(path=binary-3d.pdgm)

Load the computed diagrams.

pdlist = hc.PDList("binary-3d.pdgm")

Analysis of 1-dimensional persistence diagram

We'll examine the 1-dim persistence diagram.

pd1 = pdlist.dth_diagram(1)

pd1.histogram((-15.5, 10.5), 26).plot(colorbar={"type": "log"})

It looks like there are two birth-death pairs at (-4, 4).

pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs

[Pair(-4.0, 4.0), Pair(-4.0, 4.0)]

We have extracted two pairs.

Inverse Analysis

Let's examine what kind of holes these pairs represent.
Since we're dealing with 1-dimensional persistent homology, we should be able to identify some kind of ring structure. We'll compute this.

Here, we'll use a feature called optimal 1-cycle. This feature uses a shortest path algorithm to search for the ring with the minimal length.
It can be computed using Pair.optimal_1_cycle.

optimal_1_cycles = [pair.optimal_1_cycle() for pair in pairs]

We'll visualize them. We'll overlay the original voxel data on the optimal 1-cycles.

fig = go.Figure(data=[
    cycle.to_plotly3d() for cycle in optimal_1_cycles
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap"),
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.1, selector=dict(name="bitmap"))
fig.show()

Let's try zooming in, rotating, and adjusting the opacity and color to observe the structure. The two extracted rings surround tunnel-like structures in the voxel data.

However, the way they enclose the tunnels seems somewhat inaccurate. This is because the computation based on the shortest path is only an approximation.

Next, let's check the path that of this ring structure.
By using Optimal1CycleForBitmap.path, we can obtain the coordinates of the pixels it passes through.

optimal_1_cycles[0].path()

[(49, 4, 44),
 (49, 4, 43),
 (49, 4, 42),
 (49, 4, 41),
 (49, 4, 40),
 (49, 4, 39),
 (49, 4, 38),
 (48, 4, 38),
 (48, 4, 37),
 (48, 4, 36),
 (47, 4, 36),
 (47, 4, 35),
 (47, 4, 34),
 (46, 4, 34),
 (45, 4, 34),
 (44, 4, 34),
 (43, 4, 34),
 (43, 4, 33),
 (42, 4, 33),
 (41, 4, 33),
 (40, 4, 33),
 (39, 4, 33),
 (38, 4, 33),
 (37, 4, 33),
 (36, 4, 33),
 (36, 5, 33),
 (35, 5, 33),
 (34, 5, 33),
 (34, 6, 33),
 (34, 7, 33),
 (34, 8, 33),
 (34, 9, 33),
 (34, 10, 33),
 (34, 11, 33),
 (34, 11, 34),
 (34, 11, 35),
 (34, 11, 36),
 (34, 11, 37),
 (34, 11, 38),
 (34, 11, 39),
 (34, 11, 40),
 (34, 11, 41),
 (34, 11, 42),
 (34, 11, 43),
 (34, 11, 44),
 (34, 11, 45),
 (35, 11, 45),
 (35, 11, 46),
 (36, 11, 46),
 (37, 11, 46),
 (37, 11, 47),
 (37, 11, 48),
 (38, 11, 48),
 (39, 11, 48),
 (40, 11, 48),
 (40, 10, 48),
 (41, 10, 48),
 (42, 10, 48),
 (43, 10, 48),
 (43, 9, 48),
 (44, 9, 48),
 (45, 9, 48),
 (45, 8, 48),
 (45, 8, 47),
 (46, 8, 47),
 (46, 7, 47),
 (46, 7, 46),
 (47, 7, 46),
 (47, 7, 45),
 (48, 7, 45),
 (48, 7, 44),
 (48, 6, 44),
 (48, 5, 44),
 (49, 5, 44),
 (49, 4, 44)]

optimal_1_cycles[1].path()

[(16, 42, 45),
 (16, 41, 45),
 (16, 41, 46),
 (15, 41, 46),
 (15, 40, 46),
 (15, 40, 47),
 (14, 40, 47),
 (14, 39, 47),
 (14, 39, 48),
 (13, 39, 48),
 (13, 39, 49),
 (12, 39, 49),
 (11, 39, 49),
 (11, 38, 49),
 (11, 37, 49),
 (10, 37, 49),
 (10, 36, 49),
 (10, 35, 49),
 (10, 34, 49),
 (9, 34, 49),
 (9, 33, 49),
 (9, 32, 49),
 (8, 32, 49),
 (8, 31, 49),
 (8, 30, 49),
 (8, 29, 49),
 (8, 29, 48),
 (8, 29, 47),
 (8, 29, 46),
 (8, 29, 45),
 (8, 29, 44),
 (8, 29, 43),
 (8, 29, 42),
 (8, 29, 41),
 (9, 29, 41),
 (9, 29, 40),
 (9, 29, 39),
 (10, 29, 39),
 (10, 29, 38),
 (10, 29, 37),
 (10, 29, 36),
 (10, 29, 35),
 (10, 30, 35),
 (10, 31, 35),
 (10, 32, 35),
 (10, 33, 35),
 (11, 33, 35),
 (11, 34, 35),
 (12, 34, 35),
 (13, 34, 35),
 (13, 35, 35),
 (14, 35, 35),
 (14, 36, 35),
 (14, 37, 35),
 (14, 38, 35),
 (14, 39, 35),
 (14, 40, 35),
 (14, 41, 35),
 (14, 41, 36),
 (15, 41, 36),
 (15, 41, 37),
 (15, 42, 37),
 (15, 42, 38),
 (16, 42, 38),
 (16, 42, 39),
 (16, 42, 40),
 (17, 42, 40),
 (18, 42, 40),
 (18, 42, 41),
 (18, 42, 42),
 (18, 42, 43),
 (18, 42, 44),
 (18, 42, 45),
 (17, 42, 45),
 (16, 42, 45)]

Analysis of the 0-Dimensional Persistence Diagram

Next, let's analyze the 0-dimensional persistence diagram.

We'll plot the 0-dimensional persistence diagram.

pd0 = pdlist.dth_diagram(0)
pd0.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})

Let's take a closer look at the two pairs: (-16.0, -10.0) and (-15.0, -10.0).

For inverse analysis of 0-dimensional and 2-dimensional persistence diagrams of voxel data, we use BitmapPHTrees class (since it cannot be used for 1-dimensional data, we used optimal_1_cycle as explained above instead).

BitmapPHTrees.for_bitmap_levelset computes the information needed for inverse analysis with an interface similar to PDList.from_bitmap_levelset.

The name "PHTrees" indicates that the computed structures are trees. For n-dimensional bitmap data, it simultaneously computes two types of information: 0-dimensional and (n-1)-dimensional, resulting in two tree structures.
To extract these tree structures, we use the method PDList.bitmap_phtrees.

hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(pict, signed=True), 
                                     save_to="binary-3d-tree.pdgm")

PDList(path=binary-3d-tree.pdgm)

Load the tree structures.

pdlist_with_tree = hc.PDList("binary-3d-tree.pdgm")

Use bitmap_phtrees(0) to retrieve the tree structures from 0-dim persistent homology.

phtree_0 = pdlist_with_tree.bitmap_phtrees(0)

To obtain information on the two pairs (-16.0, -10.0) and (-15.0, -10.0), use pair_nodes_in_rectangle. This extracts all pairs (corresponding to the tree structure nodes) within the specified rectangular range.

nodes_0 = phtree_0.pair_nodes_in_rectangle(-16, -15, -10, -10)
nodes_0

[BitmapPHTrees.Node(-15.0, -10.0), BitmapPHTrees.Node(-16.0, -10.0)]

Each data is an instance of the BitmapPHTrees.Node class. We will visualize them.

fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_0
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

Observe the visualized information carefully. Changing the opacity will be helpful.

Since the data come from 0th persistent homology, the two nodes correspond to the areas inside the input voxel.

In fact, these two pairs are parent-child relationships in the tree structure, and one is a subset of the other.

The red areas may appear to have an unnatural shape.

This can be improved by using stable volume instead. For binary images, specifying stable_volume(k) will reduce the area by $k$ pixels and extract only the main part (the process is different for grayscale images). Based on experience, setting $k=1$ often works well.

fig = go.Figure(data=[
    n.stable_volume(1).to_plotly3d() for n in nodes_0
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

Each block will look more natural.

This also separates the two areas. Since stable volume extracts noise-resistant parts, this means that the parent-child relationship between the two pairs is easily disrupted by noise.

When using inverse analysis in 0-dimensional or 2-dimensional voxel data, we recommend setting it to stable_volume(1) as shown above.

Analysis of 2-dimensional persistence diagrams and inverse analysis

Finally, let's analyze the 2D persistence diagram.

Plot the 2D persistence diagram.

pd2 = pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})

(4, 6)というペアを調べます。 bitmap_ph_trees(2) として2次元パーシステンス図上の木構造を取り出します。

phtree_2 = pdlist_with_tree.bitmap_phtrees(2)

pair_nodes_in_rectangle extracts all nodes that correpond to the pairs at (4, 6).

nodes_2 = phtree_2.pair_nodes_in_rectangle(4, 4, 6, 6)
nodes_2

[BitmapPHTrees.Node(4.0, 6.0), BitmapPHTrees.Node(4.0, 6.0)]

There are two pairs of coordinates that are the same. These will be visualized in the following cell. You can see that they correspond to the centers of the voids.

fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_2
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

You can also calculate a more compact representation using stable volume.

fig = go.Figure(data=[
    n.stable_volume(1).to_plotly3d() for n in nodes_2
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

You can use the volume method to obtain the coordinates of voxels contained in this volume.

nodes_2[0].volume()

[[27, 24, 25],
 [26, 24, 25],
 [27, 23, 23],
 [26, 24, 23],
 [27, 24, 23],
 [26, 23, 23],
 [26, 23, 24],
 [27, 24, 24],
 [26, 24, 24],
 [27, 23, 24],
 [27, 23, 25],
 [26, 24, 22],
 [27, 23, 26],
 [26, 24, 26],
 [27, 24, 26],
 [28, 25, 25],
 [28, 24, 25],
 [25, 25, 25],
 [26, 25, 25],
 [27, 25, 25],
 [28, 23, 25],
 [25, 23, 25],
 [26, 23, 25],
 [25, 24, 25],
 [27, 22, 25],
 [25, 22, 25],
 [26, 22, 25],
 [26, 25, 22],
 [27, 25, 22],
 [26, 22, 23],
 [25, 22, 23],
 [27, 22, 23],
 [25, 23, 23],
 [28, 24, 23],
 [25, 24, 23],
 [28, 23, 23],
 [28, 25, 23],
 [27, 25, 23],
 [26, 25, 23],
 [25, 25, 23],
 [27, 24, 21],
 [26, 24, 21],
 [27, 23, 21],
 [26, 23, 21],
 [27, 23, 22],
 [25, 24, 22],
 [26, 23, 22],
 [27, 24, 22],
 [28, 25, 24],
 [27, 25, 24],
 [26, 25, 24],
 [28, 23, 24],
 [25, 25, 24],
 [28, 24, 24],
 [25, 24, 24],
 [28, 22, 24],
 [25, 23, 24],
 [27, 22, 24],
 [26, 22, 24],
 [25, 22, 24],
 [29, 24, 25],
 [29, 25, 25],
 [24, 25, 25],
 [25, 26, 25],
 [26, 26, 25],
 [27, 26, 25],
 [28, 26, 25],
 [26, 24, 20],
 [26, 23, 20],
 [27, 23, 20],
 [27, 24, 20],
 [27, 22, 21],
 [28, 22, 21],
 [26, 22, 21],
 [25, 22, 21],
 [24, 25, 24],
 [29, 24, 24],
 [29, 25, 24],
 [25, 26, 24],
 [26, 26, 24],
 [27, 26, 24],
 [28, 26, 24],
 [28, 23, 21],
 [28, 25, 21],
 [25, 23, 21],
 [28, 24, 21],
 [25, 24, 21],
 [25, 25, 21],
 [26, 25, 21],
 [27, 25, 21],
 [29, 22, 24],
 [24, 23, 24],
 [27, 21, 24],
 [26, 21, 24],
 [25, 21, 24],
 [24, 22, 24],
 [28, 21, 24],
 [24, 24, 24],
 [26, 22, 19],
 [29, 23, 24],
 [26, 23, 19]]

This tutorial ends here.