尤拉數

此範例說明了 2D 和 3D 物件中尤拉數 [1] 的計算。

對於 2D 物件，尤拉數是物件的數量減去孔洞的數量。請注意，如果物件考量 8 個連通像素的鄰域（2 連通性），則這相當於為互補集合（孔洞、背景）考量 4 個連通像素的鄰域（1 連通性），反之亦然。也可以使用 skimage.measure.label() 來計算物件的數量，並從兩個數之間的差推斷出孔洞的數量。

對於 3D 物件，尤拉數是透過物件的數量加上孔洞的數量，減去隧道或迴路的數量而獲得。如果對物件使用 3 連通性（將周圍 26 個體素視為其鄰域），則這對應於對互補集合（孔洞、背景）使用 1 連通性，即對於給定的體素僅考量 6 個鄰居。這裡用藍色透明表面表示體素。內部孔隙度以紅色表示。

[1]

https://en.wikipedia.org/wiki/Euler_characteristic

from skimage.measure import euler_number, label
import matplotlib.pyplot as plt
import numpy as np


# Sample image.
SAMPLE = np.array(
    [
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
        [1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0],
        [0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
        [0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    ]
)
SAMPLE = np.pad(SAMPLE, 1, mode='constant')

fig, ax = plt.subplots()
ax.imshow(SAMPLE, cmap=plt.cm.gray)
ax.axis('off')
e4 = euler_number(SAMPLE, connectivity=1)
object_nb_4 = label(SAMPLE, connectivity=1).max()
holes_nb_4 = object_nb_4 - e4
e8 = euler_number(SAMPLE, connectivity=2)
object_nb_8 = label(SAMPLE, connectivity=2).max()
holes_nb_8 = object_nb_8 - e8
ax.set_title(
    f'Euler number for N4: {e4} ({object_nb_4} objects, {holes_nb_4} '
    f'holes), \n for N8: {e8} ({object_nb_8} objects, '
    f'{holes_nb_8} holes)'
)
plt.show()

3D 物件

在此範例中，會產生一個 3D 立方體，然後加入孔洞和隧道。尤拉數會以 6 和 26 個鄰域配置進行評估。此程式碼的靈感來自 https://matplotlib.dev.org.tw/devdocs/gallery/mplot3d/voxels_numpy_logo.html

def make_ax(grid=False):
    ax = plt.figure().add_subplot(projection='3d')
    ax.grid(grid)
    ax.set_axis_off()
    return ax


def explode(data):
    """visualization to separate voxels

    Data voxels are separated by 0-valued ones so that they appear
    separated in the matplotlib figure.
    """
    size = np.array(data.shape) * 2
    data_e = np.zeros(size - 1, dtype=data.dtype)
    data_e[::2, ::2, ::2] = data
    return data_e


# shrink the gaps between voxels


def expand_coordinates(indices):
    """
    This collapses together pairs of indices, so that
    the gaps in the volume array will have a zero width.
    """
    x, y, z = indices
    x[1::2, :, :] += 1
    y[:, 1::2, :] += 1
    z[:, :, 1::2] += 1
    return x, y, z


def display_voxels(volume):
    """
    volume: (N,M,P) array
            Represents a binary set of pixels: objects are marked with 1,
            complementary (porosities) with 0.

    The voxels are actually represented with blue transparent surfaces.
    Inner porosities are represented in red.
    """

    # define colors
    red = '#ff0000ff'
    blue = '#1f77b410'

    # upscale the above voxel image, leaving gaps
    filled = explode(np.ones(volume.shape))

    fcolors = explode(np.where(volume, blue, red))

    # Shrink the gaps
    x, y, z = expand_coordinates(np.indices(np.array(filled.shape) + 1))

    # Define 3D figure and place voxels
    ax = make_ax()
    ax.voxels(x, y, z, filled, facecolors=fcolors)
    # Compute Euler number in 6 and 26 neighborhood configuration, that
    # correspond to 1 and 3 connectivity, respectively
    e26 = euler_number(volume, connectivity=3)
    e6 = euler_number(volume, connectivity=1)
    plt.title(f'Euler number for N26: {e26}, for N6: {e6}')
    plt.show()


# Define a volume of 7x7x7 voxels
n = 7
cube = np.ones((n, n, n), dtype=bool)
# Add a tunnel
c = int(n / 2)
cube[c, :, c] = False
# Add a new hole
cube[int(3 * n / 4), c - 1, c - 1] = False
# Add a hole in neighborhood of previous one
cube[int(3 * n / 4), c, c] = False
# Add a second tunnel
cube[:, c, int(3 * n / 4)] = False
display_voxels(cube)