Surface area computation

This notebook showcases various ways to compute specific surface areas and interfacial areas of labelled phases based on dummy and real-world examples.

import matplotlib.pyplot as plt
import numpy as np
import taufactor.metrics as tau
import tifffile
from scipy.ndimage import rotate

Specific surface area validation

• cube - perfectly grid-aligned

• rotated cube (45° around z-axis)

• multiple cubes

• two {111}-planes in cube

• sphere

Let’s start by creating the structures…

nx, ny, nz = [20, 20, 20]

cube = np.zeros((nx, ny, nz))
cube[5:-5,5:-5,5:-5] = 1

cube_rot = rotate(
    cube, angle=45.0, axes=(1, 0),  # rotate in (y, x) -> about z
    reshape=False,
    order=0,
    mode='constant',
    cval=0
).astype(np.uint8)

cubes = np.zeros((nx, ny, nz))
cubes[0:10, 0:10, 0:5] = 1
cubes[5:-5, 5:-5, 5:-5] = 2
cubes[0:10, 0:10, 15:] = 1
cubes[10:, 5:, 15:] = 3

diagonal = np.zeros((nx, ny, nz))
x, y, z = np.ogrid[:nx, :ny, :nz]
plane = (x+0.5)/nx + (y+0.5)/ny + (z+0.5)/nz
mask = (plane <= 1.05) | (plane > 2.0)
diagonal[mask] = 1

sphere = np.zeros((nx, ny, nz))
radius = np.min([nx,ny,nz])*0.5-3
x, y, z = np.ogrid[:nx, :ny, :nz]
distance_squared = (x - nx/2 + 0.5)**2 + (y - ny/2 + 0.5)**2 + (z - nz/2 + 0.5)**2
mask = distance_squared <= radius**2
sphere[mask] = 1

structures = {'cube': cube,
              'cube_rot': cube_rot,
              'multicube': cubes,
              'diagonal': diagonal,
              'sphere': sphere}

a_theo = {'cube': (2*(nx-10)*(ny-10)+2*(nx-10)*(nz-10)+2*(ny-10)*(nz-10))/(nx*ny*nz),
          'cube_rot': (2*(nx-10)*(ny-10)+2*(nx-10)*(nz-10)+2*(ny-10)*(nz-10))/(nx*ny*nz),
          'multicube': (2*(nx-10)*(ny-10)+(nx-10)*(nz-10)+(ny-10)*(nz-10))/(nx*ny*nz),
          'diagonal': np.sqrt(3)*nx**2/(nx*ny*nz),
          'sphere': 4*np.pi*radius**2/(nx*ny*nz)}

labels = {'gray': 1}

…and visualize the structures.

fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(16, 4), dpi=100)

i = 0
for key, value in structures.items():
    axes[i].remove()
    axes[i] = fig.add_subplot(1, 5, i + 1, projection='3d')

    axes[i].voxels(value==1, facecolors='gray', edgecolor='black', alpha=1.0)
    axes[i].voxels(value==2, facecolors='blue', edgecolor='black', alpha=1.0)
    axes[i].voxels(value==3, facecolors='green', edgecolor='black', alpha=1.0)
    axes[i].set_aspect('equal')
    axes[i].set_title(key)
    axes[i].set_xlim(0, nx)
    axes[i].set_ylim(0, ny)
    axes[i].set_zlim(0, nz)
    i = i + 1
plt.tight_layout()
plt.show()

The gray phase is labelled with 1 and will be used for the comparison of various surface areas computations. Note that already from the volume fractions we can see grid anisotropy i.e. the rotated cube now has 12% more voxels.

for key, struc in structures.items():
    print("Volume fraction of "+key+f": {tau.volume_fraction(struc, labels)['gray']:.4f}")

Volume fraction of cube: 0.1250
Volume fraction of cube_rot: 0.1400
Volume fraction of multicube: 0.1250
Volume fraction of diagonal: 0.3587
Volume fraction of sphere: 0.1840

a_faces = {}
a_marching = {}
a_gradient = {}

for key, struc in structures.items():
    a_faces[key] = tau.specific_surface_area(struc,
                                             method='face_counting',
                                             phases=labels)['gray']
    a_gradient[key] = tau.specific_surface_area(struc,
                                                method='gradient',
                                                phases=labels,
                                                smoothing=False)['gray']
    a_marching[key] = tau.specific_surface_area(struc,
                                                method='marching_cubes',
                                                phases=labels,
                                                smoothing=False,
                                                device='cpu')['gray']

print("Spec. surf. area:     faces         \t|      gradient      \t|   marching")
for key, struc in structures.items():
    print(key+f":   \t{a_faces[key]:.5f} ({((a_faces[key]-a_theo[key])/a_theo[key]*100):.2f} %)\t|"
        f"  {a_gradient[key]:.5f} ({(a_gradient[key]-a_theo[key])/a_theo[key]*100:.2f} %)\t|"
        f"  {a_marching[key]:.5f} ({(a_marching[key]-a_theo[key])/a_theo[key]*100:.2f} %)")

Spec. surf. area:     faces             |      gradient         |   marching
cube:           0.07500 (0.00 %)        |  0.07085 (-5.53 %)    |  0.07051 (-5.98 %)
cube_rot:       0.09800 (30.67 %)       |  0.07657 (2.10 %)     |  0.07624 (1.65 %)
multicube:      0.05000 (0.00 %)        |  0.04823 (-3.54 %)    |  0.04219 (-15.61 %)
diagonal:       0.14963 (72.77 %)       |  0.09402 (8.57 %)     |  0.07813 (-9.78 %)
sphere:         0.11700 (52.01 %)       |  0.08378 (8.85 %)     |  0.08344 (8.41 %)

From the comparison of surface area computations we can see that face counting is exact for all cases where we have grid aligned structures (cube and multicube). The specific surface area of a sphere is overestimated by \(50\%\) even in the limit of high voxel resolution as discussed by Daubner et. al. and, similarly, the diagonal planes are overestimated by \(\sqrt{3}-1 = 73.2\%\).

The marching_cubes algorithm and the gradient method introduce some errors for grid aligned cases (due to edge smoothing) but generally perform much better for curved surfaces as often encountered in porous microstructures. Especially in the context of battery electrodes, the spherical example is the most relevant comparison. Both the marching_cubes and the gradient method can be combined with a smoothing operation (convolution with a 3x3x3 Gaussian kernel) which further increases the accuracy for curved surfaces but leads to increased edge smoothing for all cuboid structures. Note that smoothing=True is the default option for the specific_surface_area function.

for key, value in structures.items():
    a_gradient[key] = tau.specific_surface_area(value, method='gradient', phases=labels)['gray']
    a_marching[key] = tau.specific_surface_area(value, method='marching_cubes', phases=labels, device='cpu')['gray']

print("Spec. surf. area:     gradient      \t|   marching")
for key, struc in structures.items():
    print(key+f":     \t{a_gradient[key]:.5f} ({(a_gradient[key]-a_theo[key])/a_theo[key]*100:.2f} %)\t|"
        f"  {a_marching[key]:.5f} ({(a_marching[key]-a_theo[key])/a_theo[key]*100:.2f} %)")

Spec. surf. area:     gradient          |   marching
cube:           0.06547 (-12.71 %)      |  0.06368 (-15.10 %)
cube_rot:       0.07063 (-5.83 %)       |  0.06941 (-7.45 %)
multicube:      0.04587 (-8.26 %)       |  0.03958 (-20.83 %)
diagonal:       0.08502 (-1.83 %)       |  0.07626 (-11.94 %)
sphere:         0.07770 (0.95 %)        |  0.07634 (-0.82 %)

Interfacial area validation

Next, let’s look at pairwise interfacial areas. For now, only face_counting is implemented.

In the case of a binary structure, there will only be one interfacial area which must be equivalent to to respective specific surface areas. This can be easily shown:

print(tau.specific_surface_area(sphere, method='face_counting'))
print(tau.interfacial_areas(sphere, normalize=True))

{'0': 0.117, '1': 0.117}
{(0, 1): 0.117}

The multi-label cube is a good validation example for this function. Note that there is a normalize = True option to compute the specific interfacial areas in [1/m] while normalize = False returns the toal amount of counted faces.

tau.interfacial_areas(cubes, normalize = False)

{(0, 1): 325.0,
 (0, 2): 500.0,
 (0, 3): 200.0,
 (1, 2): 50.0,
 (1, 3): 25.0,
 (2, 3): 50.0}

These results are correct which can be easily confirmed by looking at the multi-cube structure in the figure above.

Battery electrode example

structure = tifffile.imread('electrode.tiff')
print("Stack shape:", structure.shape)
labels = {"pore":0, "NMC":85, "CBD":170}

Stack shape: (256, 256, 256)

Let’s track the wall time of each method

import time
timer = []
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='face_counting', device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='gradient', smoothing=False, device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='gradient', smoothing=True, device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='marching_cubes', smoothing=False, device='cpu'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='marching_cubes', smoothing=True, device='cpu'))
timer.append(time.time())

print(f"Face_counting: {timer[1]-timer[0]:.4f} s\n"
      f"Gradient w/o smoothing: {timer[2]-timer[1]:.4f} s\n"
      f"Gradient w smoothing: {timer[3]-timer[2]:.4f} s\n"
      f"Marching w/o smoothing: {timer[4]-timer[3]:.4f} s\n"
      f"Marching w smoothing: {timer[5]-timer[4]:.4f} s"
      )

{'pore': 0.30683690309524536, 'NMC': 0.14610451459884644, 'CBD': 0.34222638607025146}
{'pore': 0.19837072491645813, 'NMC': 0.10425256937742233, 'CBD': 0.2104557752609253}
{'pore': 0.1161952018737793, 'NMC': 0.09140095114707947, 'CBD': 0.109257273375988}
{'pore': 0.22268114984035492, 'NMC': 0.10418400168418884, 'CBD': 0.25399911403656006}
{'pore': 0.136002317070961, 'NMC': 0.09128409624099731, 'CBD': 0.13165849447250366}
Face_counting: 0.2399 s
Gradient w/o smoothing: 0.1088 s
Gradient w smoothing: 0.1400 s
Marching w/o smoothing: 7.5878 s
Marching w smoothing: 6.1020 s

From these results we can already see that the marching cubes method is more than an order of magnitude slower than the other methods.

Next we use the SNOW segmentation to identify all individual NMC particles

import porespy as ps

snow_labels = (ps.filters.snow_partitioning(structure == labels["NMC"])).regions
print(f"Identified {np.unique(snow_labels).size} individual particles.")

Identified 1363 individual particles.

timer = []
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='face_counting', device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='gradient', smoothing=False, device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='gradient', smoothing=True, device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='marching_cubes', smoothing=False, device='cpu')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='marching_cubes', smoothing=True, device='cpu')
timer.append(time.time())

print(f"Face_counting: {timer[1]-timer[0]:.4f} s\n"
      f"Gradient w/o smoothing: {timer[2]-timer[1]:.4f} s\n"
      f"Gradient w smoothing: {timer[3]-timer[2]:.4f} s\n"
      f"Marching w/o smoothing: {timer[4]-timer[3]:.4f} s\n"
      f"Marching w smoothing: {timer[5]-timer[4]:.4f} s"
      )

Face_counting: 0.4330 s
Gradient w/o smoothing: 3.3783 s
Gradient w smoothing: 4.0257 s
Marching w/o smoothing: 74.9045 s
Marching w smoothing: 76.2564 s

interfaces = tau.interfacial_areas(snow_labels, normalize=False)
count = 0
av_interface_area = 0
for pair, val in interfaces.items():
    if pair[0] != 0 and pair[1] != 0:
        count += 1
        av_interface_area += val
av_interface_area = av_interface_area / count
print(f"Average interfacial area between particles is {av_interface_area} voxel faces.")

Average interfacial area between particles is 111.51503957783642 voxel faces.