3D binary image analysis

This tutorial explains how to analyze a 2D binary image by persistent homology. The basic idea is the same as 2D binary image analysis, so please try this tutorial after you have finished that tutorial.

Here, 50x50x50 size image is analyzed. There are various ways to represent 3D images, but here we use cross-sectional images for each layer. In other words, prepare 50x50x50 images. The file is located in the 'data/' directory. Name the file in the order of the 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

Read data file

First, read this 50x50x50 image and make a 50x50x50 black and white array.

# Loading required libraries
import imageio
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Before creating a 3D array, load 0000.png and check what image it is.

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

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

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

This is a 50x50 array with values from 0 to 255. So this image is like a grayscale image. Let's look at the contents and actually display them.

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]], dtype=uint8)

plt.imshow(image_0000, "gray")

<matplotlib.image.AxesImage at 0x7ff23affa530>

Now load 50 images 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 50x50x50 image has been completed. Let's visualize this in 3D. We use pyvista. A translucent representation is realized by opacity=0.7.

import pyvista as pv  # Load pyvsta module
import homcloud.pyvistahelper as pvhelper  # HomCloud's pyvista helpers

pv.set_jupyter_backend('trame') 

pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.7)
pl.show()

Widget(value="<iframe src='http://localhost:35151/index.html?ui=P_0x7ff2b41a7730_0&reconnect=auto' style='widt…

Calculation of persistence diagram

Calculate persistence diagram. Use hc.PDList.from_bitmap_distance_function. Same as 2D.

import homcloud.interface as hc

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

You can also load the file as follows.

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

Analysis of persistence diagram (1st)

Let's look at the first persistence diagram.

pd1 = pdlist.dth_diagram(1)

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

There seems to be two birth-death pairs at (-4, 4). Let's check this.

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

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

There are certainly two.

Inverse analysis(1st)

Let's take a look at what the two pairs represent. Because you are looking at the first-order persistence homology, you should be looking at some kind of ring structure. Let's calculate it.

Here you use the optimal cycle. In point cloud analysis, we used optimal volume and volume-optimal cycle. It's similar (slightly different) Calculate with Pair.optimal_1_cycle.

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

Let's visualize this first. Displays the original voxel data superimposed.

pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.7)
for cycle in optimal_1_cycles:
    pl.add_mesh(cycle.to_pyvista_mesh(), color="red")
pl.show()

Widget(value="<iframe src='http://localhost:35151/index.html?ui=P_0x7ff1c78cf3d0_1&reconnect=auto' style='widt…

Please operate Paraview in various ways. Change the Opacity of voxel data in various ways, increase the thickness of optimal 1-cycles line, etc. Explore settings that are easy to see. As you can see from visualization, in the case of voxel data, the tunnel structure corresponds to the first-order birth-death pair.

Next, let's see what path this ring structure is taking. Use Optimal1CycleForBitmap.path. You can get the coordinates of pixel at the passing point.

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 persistence diagram(0th) and its inverse analysis

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

First, plot a 0th-order 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 birth-death pairs, (-16.0, -10.0), (-15.0, -10.0). Use a class called BitmapPHTrees for the 0th and 2nd order inverse analysis of voxel data(it cannot be used for the 1st order, so we used the optimal_1_cycle described above instead). BitmapPHTrees.for_bitmap_levelset has an interface similar to PDList.from_bitmap_levelset, calculate information for back analysis.

The name PHTrees means that what is computed here is a tree structure. BitmapPHTrees.for_bitmap_levelset simulateneously computes two forests for the information of 0th and (n-1)th PH of an n dimensional bitmap.

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 data calculated as follows.

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

Use bitmap_phtrees(0) to retrieve the 0th PH's information.

phtree_0 = pdlist_with_tree.bitmap_phtrees(0)

Use pair_nodes_in_rectangle to retrieve information on the two birth-death pairs (-16.0, -10.0), (-15.0, -10.0). This retrieves all birth-death pairs (with corresponding tree nodes) within the specified 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)]

The extracted data is an instance of the class called BitmapPHTrees.Node. Next, let's visualize this information.

pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_0[0].to_pyvista_mesh(), color="green")
pl.add_mesh(nodes_0[1].to_pyvista_mesh(), color="red", opacity=0.4)
pl.show()

Widget(value="<iframe src='http://localhost:35151/index.html?ui=P_0x7ff1c78f32e0_2&reconnect=auto' style='widt…

Since it is the 0th PH, it corresponds to the area inside the voxel. Actually, these two birth-death pairs are in a parent-child relationship on a tree structure, and one is a subset of the other.

The green area looks jagged, and something unnatural. In this example, we can solve the problem using stable volume instead. In many cases, it looks more natural if you use stable_volume(1) to shrink the volume by 1 voxel, as shown below.

pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_0[0].stable_volume(1).to_pyvista_mesh(), color="green")
pl.add_mesh(nodes_0[1].stable_volume(1).to_pyvista_mesh(), color="red", opacity=0.4)
pl.show()

Widget(value="<iframe src='http://localhost:35151/index.html?ui=P_0x7ff1c894ace0_3&reconnect=auto' style='widt…

Each clump looks better than before, and the two regions are now separated. This phenomenon means that the parent-child relationship between two birth-death pairs is easily corrupted by noise since the mechanism of stable volumes removes the noise-sensitive part of the volume.

Analysis of persistence diagram (secondary) and its inverse analysis

Finally, let's analyze the second-order persistence diagram.

Plot a second-order persistence diagram.

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

Examine the birth-death pair called (4, 6). For inverse analysis, use bitmap_phtrees(2).

phtree_2 = pdlist_with_tree.bitmap_phtrees(2)

Let's look at the birth-death pair called (4, 6). First, extract all birth-death pairs (4, 6)(nodes of the tree corresponding to) with pair_nodes_in_rectangle.

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)]

Two the following is the result of visualization.　You can see that the deeper part of the hollow area corresponds to this birth-death pair.

pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_2[0].to_pyvista_mesh())
pl.add_mesh(nodes_2[1].to_pyvista_mesh())
pl.show()

Widget(value="<iframe src='http://localhost:35151/index.html?ui=P_0x7ff1c78ce470_4&reconnect=auto' style='widt…

You can also get each coordinate of this area by using the volume method.

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]]