Minkowski Functionals#

This notebook demonstrates some methods of computing the four Minkowski functionals (in 3D). The four Minkowski functionals are:

  1. Volume

  2. Surface Area

  3. Integral Mean Curvature

  4. Euler Characteristic

In this notebook, we compare Minkowski functional calculations via 2 methods:

Import some packages#

# Some Utility functions
from pathlib import Path

# Image processsing
import skimage
from skimage.morphology import ball

# Image visualization
import pyvista as pv
pv.set_jupyter_backend('static')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

# Image quantification
import numpy as np
import porespy as ps
import trimesh
from dpm_tools.io import read_image, Image
from dpm_tools.metrics import minkowski_functionals
from dpm_tools.visualization import plot_isosurface

1. Computing the Minkowski Functionals#

For the following introductory demonstrations, we will use a single solid sphere sample as the motivating example. At the end, we will have the opportunity to try the functions on the other samples.

Generally speaking, the Minkowski functionals can be computed using voxel methods or mesh methods. We will explore both here using the Quantimpy (voxelized) and the Trimesh/Skimage (mesh) libraries.

In the following cell, we create and visualize a solid sphere of radius 100 voxels

radius = 100
side_length = 250

sphere = np.zeros([side_length, side_length, side_length], dtype=np.uint8)
id_a = side_length // 2 - radius
id_b = side_length // 2 + radius + 1
sphere[id_a:id_b, id_a:id_b, id_a:id_b] = ball(radius, dtype=np.uint8)

# DPM Tools Image object
sphere = Image(scalar=sphere)

Before we begin, we create a function to compare our package computations with the analytical solution

def analytical_solution(radius):
    volume_analytical = 4/3 * np.pi * radius**3
    surface_area_analytical = 4 * np.pi * radius**2
    mean_curvature_analytical = 4 * np.pi * radius
    euler_analytical = 1.0
    
    return volume_analytical, surface_area_analytical, mean_curvature_analytical, euler_analytical


def compute_error(true, computed_measures):
    return [np.abs((true - compute) / true)*100 for compute in computed_measures]

volume_analytical, surface_area_analytical, mean_curvature_analytical, euler_analytical = analytical_solution(radius)

1.1 Volume Calculation#

Considering a 3D body, \(X\), with a smooth boundary, \(\delta X\), the first Minkowski functional (volume) can be computed as:

\[ M_0(X) = \int_X dv \]

For a sphere, the volume can be found as:

\[ M_0 = \frac{4}{3} \pi R^3 \]

In practice, volume is the most straightforward of the Minkowski functionals set to compute.

Here, we compare a few different methods and packages for computing the volume of the pore space, split into voxelized methods and mesh methods:

  1. Voxel Counting

  2. Quantimpy (Voxelized)

  3. Trimesh (Mesh)

def compute_volume(image):
    # Simple voxel counting
    voxel_counting = np.sum(image == 1)

    # DPM Tools
    dpm_tools_volume = minkowski_functionals(image)[0]
    
    # Mesh Volume
    mesh_measures = np.abs(ps.metrics.mesh_volume(image))
    
    return voxel_counting, dpm_tools_volume, mesh_measures

volume_measures = compute_volume(sphere.scalar)
volume_errors = compute_error(volume_analytical, volume_measures)

print("Volume Calculations:")
print(f"\tAnalytical Solution: {volume_analytical:0.3f}")
print(f"\tVoxel Counting: {volume_measures[0]:0.3f} \t Rel. Error: {volume_errors[0]:0.3f}%")
print(f"\tDPM Tools: {volume_measures[1]:0.3f} \t Rel. Error: {volume_errors[1]:0.3f}%")
print(f"\tMesh: {volume_measures[2]:0.3f} \t\t Rel. Error: {volume_errors[2]:0.3f}%")

Volume Calculations:
	Analytical Solution: 4188790.205
	Voxel Counting: 4187857.000 	 Rel. Error: 0.022%
	DPM Tools: 4187857.000 	 Rel. Error: 0.022%
	Mesh: 4187199.480 		 Rel. Error: 0.038%

1.2 Surface Calculation#

The second Minkowski functional (surface area) can be computed as:

\[ M_1(X) = \int_{\delta X} ds \]

For a sphere, the surface area can be found as:

\[ M_1 = 4 \pi R^2 \]

Again, there are multiple ways to compute the surface area of the solid/pore boundary. Here, we compare:

  1. Quantimpy (voxelized), and

  2. Surface Mesh method (Skimage):

def compute_surface_area(image, smoothing=None):

    # DPM Tools
    dpm_tools_surface_area = minkowski_functionals(image)[1]
    
    # Mesh Surface Area
    verts, faces, normals, values = skimage.measure.marching_cubes(image, level = smoothing)
    mesh_surface_area = skimage.measure.mesh_surface_area(verts, faces)
    
    return dpm_tools_surface_area, mesh_surface_area

# Compute the surface area using 0.5 isosurface
surface_area_measures = compute_surface_area(sphere.scalar)
surface_area_errors = compute_error(surface_area_analytical, surface_area_measures)

# Compute the surface area using some other isosurface
_, surface_area_measures_blocky = compute_surface_area(sphere.scalar, smoothing:=0)
surface_area_errors_blocky = compute_error(surface_area_analytical, [surface_area_measures_blocky])

print("Surface Area Calculations:")
print(f"\tAnalytical Solution: {surface_area_analytical:0.3f}")
print(f"\tDPM_Tools: {surface_area_measures[0]:0.3f} \t\t Rel. Error: {surface_area_errors[0]:0.3f}%")
print(f"\tMesh (0.5 Isosurface): {surface_area_measures[1]:0.3f} \t Rel. Error: {surface_area_errors[1]:0.3f}%")
print(f"\tMesh ({smoothing} Isosurface): {surface_area_measures_blocky:0.3f} \t Rel. Error: {surface_area_errors_blocky[0]:0.3f}%")

Surface Area Calculations:
	Analytical Solution: 125663.706
	DPM_Tools: 125668.000 		 Rel. Error: 0.003%
	Mesh (0.5 Isosurface): 136694.734 	 Rel. Error: 8.778%
	Mesh (0 Isosurface): 137776.469 	 Rel. Error: 9.639%

1.3 Integral Mean Curvature Calculation#

The third Minkowski functional (integral mean curvature area) can be computed as:

\[ M_2(X) = \int_{\delta X} \frac{1}{2} \left[ \frac{1}{R_1} + \frac{1}{R_2} \right] ds \]

For a sphere, the principal radii of curvature are the same everywhere (i.e. \(R_1 = R_2\)). Therefore, the analytical integral mean curvature can be found as:

\[ M_2 = 4 \pi R \]

3D curvature measurements are not trivial to compute. Though there are many methods to compute the curvature, we will stick with the built in functions available in Quantimpy and the Trimesh library.

  1. Quantimpy (voxelized), and

  2. Surface Mesh method (Implemented with Trimesh):

def compute_mean_curvature(image):
    # Quantimpy
    dpm_tools_mean_curvature = minkowski_functionals(image)[2]
    
    # Mesh Surface Area
    trimesh_sphere = trimesh.creation.icosphere(radius=radius)
    mesh_mean_curvature = trimesh_sphere.integral_mean_curvature
    
    return dpm_tools_mean_curvature, mesh_mean_curvature

mean_curvature_measures = compute_mean_curvature(sphere.scalar)
mean_curvature_errors = compute_error(mean_curvature_analytical, mean_curvature_measures)

print("Mean Calculations:")
print(f"\tAnalytical Solution: {mean_curvature_analytical:0.3f}")
print(f"\tDPM Tools: {mean_curvature_measures[0]:0.3f} \t Rel. Error: {mean_curvature_errors[0]:0.3f}%")
print(f"\tMesh: {mean_curvature_measures[1]:0.3f} \t\t Rel. Error: {mean_curvature_errors[1]:0.3f}%")

Mean Calculations:
	Analytical Solution: 1256.637
	DPM Tools: 1262.920 	 Rel. Error: 0.500%
	Mesh: 1254.639 		 Rel. Error: 0.159%

1.4 Euler Characteristic#

The fourth Minkowski functional (Gaussian curvature) can be computed as:

\[ M_3(X) = \int_{\delta X} \left[ \frac{1}{R_1R_2} \right] ds \]

Because the principal radii of curvature are the same everywhere for a sphere (i.e. \(R_1 = R_2\)), the analytical total curvature can be found as:

\[ M_3 = 4 \pi\]

The Gauss-Bonnet theorem links the Gaussian Curvature to the Euler characteristic (or Euler number) by:

\[ 4\pi\chi = M_3(X).\]

So, for a solid ball, the Euler number (\(\chi\)) \(= 1\)

Here, we compute the Euler number using:

  1. Quantimpy, and

  2. Skimage

def compute_euler_number(image):
    # Quantimpy
    dpm_tools_euler_number = minkowski_functionals(image)[3]
    
    # Mesh Surface Area
    skimage_euler_number = skimage.measure.euler_number(image, connectivity=3)
    
    return dpm_tools_euler_number, skimage_euler_number

euler_number_measures = compute_euler_number(sphere.scalar)
euler_number_errors = compute_error(euler_analytical, euler_number_measures)

print("Euler Number Calculations:")
print(f"\tAnalytical Solution: {euler_analytical:0.3f}")
print(f"\tDPM Tools: {euler_number_measures[0]:0.3f} \t Rel. Error: {euler_number_errors[0]:0.3f}%")
print(f"\tSkimage: {euler_number_measures[1]:0.3f} \t\t Rel. Error: {euler_number_errors[1]:0.3f}%")

Euler Number Calculations:
	Analytical Solution: 1.000
	DPM Tools: 1.000 	 Rel. Error: 0.000%
	Skimage: 1.000 		 Rel. Error: 0.000%

2. Minkowski Functionals in Rock Samples#

In the following cell, we compute the Minkowski functionals for the four samples we previously examined. For sake of simplicity, we will only use the Quantimpy library here.

2.1 Loading in our data#

We load in TIFF images from the data directory. This directory contains segmented x-ray microtomography data from the Network Generation Comparison Forum and is available on the Digital Rocks Portal (https://www.digitalrocksportal.org/projects/16).

For the purposes of this workshop, we have preprocessed the data to be in TIFF file format with solid and pore labeled as 0 and 1, respectively.

datapath = Path('../../../_static')
beadpack = read_image(datapath / 'beadpack.tif')
gambier = read_image(datapath / 'mtgambier.tif')

fig = plt.figure(figsize=(12, 5))
gs = GridSpec(1, 3, width_ratios=[1, 1, 0.05], wspace=0.3)

# Create subplots
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
cax = fig.add_subplot(gs[0, 2])


im1 = ax1.imshow(beadpack[0,:,:], cmap='gray')
im2 = ax2.imshow(gambier[0,:,:], cmap='gray')

fig.colorbar(im1, cax=cax)

ax1.set_title('Glass Bead Pack',fontsize=14)
ax2.set_title('Mt. Gambier Limestone',fontsize=14)

plt.show()
../../../_images/2df0287f23ca1b4b320b68ce615e72c3567916ee0fc24963f75f3c28fc910493.png 
def compute_sample_mf(sample):
    
    sample_dict = {'Mt. Gambier Limestone': gambier, 
                   'Beadpack': beadpack}
    
    fig, ax = plt.subplots(nrows=1, ncols=1,figsize=(8,8))
    
    plt.imshow(sample_dict[sample][0], cmap='gray')
    ax.set_title(sample,fontsize=14)
    plt.colorbar(fraction=0.046, pad=0.04)
    mfs = minkowski_functionals(sample_dict[sample])

    print('Minkowski Functionals:')
    print(f'\tVolume: {mfs[0]:0.3f}')
    print(f'\tSurface Area: {mfs[1]:0.3f}')
    print(f'\tIntegral Mean Curvature: {mfs[2]:0.3f}')
    print(f'\tEuler Number: {mfs[3]:0.3f}')

    

compute_sample_mf('Beadpack')

Minkowski Functionals:
	Volume: 6356074.000
	Surface Area: 1432277.333
	Integral Mean Curvature: -19064.231
	Euler Number: -1455.500
../../../_images/fa77fbef638b3d31b2e4fda6efa7602f4d6424a2a25866af8477cf2ac97a9a04.png 
compute_sample_mf('Mt. Gambier Limestone')

Minkowski Functionals:
	Volume: 7314939.000
	Surface Area: 2163465.333
	Integral Mean Curvature: 118698.795
	Euler Number: -2971.500
../../../_images/58310195ef2407e1065a16d9403e06225cb3c030b0ab05ee6ff0236cd237771b.png