Structure Tensor Fiber Analysis Demo

Author: Niels Jeppesen (niejep@dtu.dk)

In this notebook we use Structure Tensor (ST) for orientation analysis of glass fiber composites. We use orientation information to both segment fibers and extract statistical information about the orientation of the fibers.

The structure-tensor package we will be using here is a 2D and 3D strcture tensor package for Python implemented by Vedrana A. Dahl and Niels Jeppesen.

Installation

To run this notebook there are some prerequisites that must be installed in our Python environment. We generally recommend creating a new Python environment to run the notebooks in, for instance using conda create -n <env_name> python=3.7. This way we avoid conflicts with packages in your existing environment.

Option 1

Install dependencies (for this notebook) using pip install numpy scipy scikit-image matplotlib nibabel structure-tensor. To start the notebook we will of course also need jupyter.

Option 2

Run pip install -r requirements.txt using the included requirements.txt file. This will install specific versions of all packages needed to run this notebook, overwriting any current versions of those packages in our environment. This also includes packages used in our other notebook.

Now, let's go ahead and import the required modules. The structure_tensor_workers.py file should be in the notebook folder.

In [1]:
import os

import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import nibabel as nib
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
from scipy import ndimage
from scipy.spatial.transform import Rotation as R
from skimage import filters, transform
from structure_tensor import eig_special_3d, structure_tensor_3d

from structure_tensor_workers import calculate_angles

Data

In this notebook, we will be working with a 3D x-ray CT scan of a non-crimp glass fiber material. The original TXM data file was cropped and converted to NIfTI (.nii) to reduce size and ease data loading.

First we need to set the file name and path. We will be using the following folder structure:

  • ../notebooks contains the notebooks.
  • ../originals should contain the original data files.
  • ../tmp should contain any temperary files or output we create.
In [2]:
# Set file name and path. You can change these to load your own NIfTI data.
file_name = 'HF401TT-13_FoV16.5_Stitch.txm.nii'
file_path = os.path.join('../originals/', file_name)
print('File path:', file_path)

File path: ../originals/HF401TT-13_FoV16.5_Stitch.txm.nii

Next step is to load the file using the nibabel package.

In [3]:
# Load data.
nii_file = nib.load(file_path)

# Read meta data.
data_shape = nii_file.shape
data_type = nii_file.get_data_dtype()
voxel_size = nii_file.affine[0, 0]

# Print meta data.
print('Shape:', data_shape)
print('Data type:', data_type)
print('Voxel size:', voxel_size, 'μm')

Shape: (596, 1922, 4228)
Data type: uint16
Voxel size: 8.078731536865234 μm

The object returned by nibabel.load allows us to access both data and metadata, but will not actually load data into memory. Instead we can access data using the dataobj property. Altertively the get_data function can be used to return data as an standard in-memory NumPy array.

Show data

After loading the data, it's always a good idea to have a look at it.

First we'll create a helper function for adding a simple scale bar to our figures.

In [4]:
def add_scalebar(ax, length=1, scale=1, unit='mm', text=None, ticks=False):
    """Add a scale bar to the plot."""
    if not ticks:
        ax.set_xticks([])
        ax.set_yticks([])
    
    text = f'{length} {unit}' if text is None else text
    
    fontprops = fm.FontProperties(size=18)
    scalebar = AnchoredSizeBar(ax.transData,
                               length / scale, text, 'lower left', 
                               pad=0.3,
                               color='white',
                               frameon=False,
                               size_vertical=3,
                               fontproperties=fontprops)
    ax.add_artist(scalebar)

Let's have a look at the center slice of all three axes.

In [5]:
fig, ax = plt.subplots(1, 1, figsize=(20, 20))
i = data_shape[2] // 2
img = nii_file.dataobj[..., i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('YZ-slice')
ax.set_xlabel('Y')
ax.set_ylabel('Z')
add_scalebar(ax, length=2, scale=voxel_size / 1000, ticks=True)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(20, 20))
i = data_shape[1] // 2
img = nii_file.dataobj[:, i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('XZ-slice')
ax.set_xlabel('X')
ax.set_ylabel('Z')
add_scalebar(ax, length=2, scale=voxel_size / 1000, ticks=True)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(20, 20))
i = data_shape[0] // 2
img = nii_file.dataobj[i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('XY-slice')
ax.set_xlabel('X')
ax.set_ylabel('Y')
add_scalebar(ax, length=2, scale=voxel_size / 1000, ticks=True)
plt.show()

Cut out data

To keep things simple we will start by cutting out a small volume to work with in this notebook. This way we don't have to worry too much about performance and memory issues. To calculate ST for larger data sets, we highly recommend splitting them into blocks during calculation. We will show how this can be done in the end of this notebook.

The data we will be working with will be stored in the data variable as a uint16 NumPy array.

In [6]:
# Cut out a small volume and put it into a NumPy array.
data = np.array(nii_file.dataobj[:, 250:500, 2000:2500])

For our sample, we want fibers to span along the X-axis, and the Z-axis to start at the "front" of the material, so we rotate the data. This is specific to the sample, so you may want to change or remove the cell below. If you change the rotation, you will also need to change the axis labels in the figures for them to be correct.

In [7]:
# Rotate to match conversion. This is specific to the data.
# If it is not relevant for you, just uncomment the line below.
data = np.rot90(data, k=3, axes=(0, 2))
In [8]:
# Print the shape of the cut-out.
data_shape = data.shape
print('Cut-out shape:', data_shape)

Cut-out shape: (500, 250, 596)

Now that we have a smaller volume to work with, let's start by once again having a look at it.

In [9]:
fig, ax = plt.subplots(1, 1, figsize=(20, 10))
i = data_shape[0] // 2
img = data[i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('YZ-slice')
ax.set_xlabel('Z')
ax.set_ylabel('Y')
add_scalebar(ax, scale=voxel_size / 1000, ticks=True)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(20, 10))
i = data_shape[1] // 2
img = data[:, i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('XZ-slice')
ax.set_xlabel('Z')
ax.set_ylabel('X')
add_scalebar(ax, scale=voxel_size / 1000, ticks=True)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(20, 10))
i = data_shape[2] // 2
img = data[..., i]
ax.imshow(img.squeeze(), cmap='gray')
ax.set_title('XY-slice')
ax.set_xlabel('Y')
ax.set_ylabel('X')
add_scalebar(ax, scale=voxel_size / 1000, ticks=True)
plt.show()

Rotate data

In cases where we want to rotate the data arbitrarily, the scipy.ndimage.affine_transform function can be used. To do so, we first need to define a rotation matrix, for instance using Euler angles as shown below. As we don't need to rotate the our data in this demo, we've commented out the code.

In [10]:
# Define rotation matrix.
# rot_matrix = R.from_euler('xyz', [0.39, 0.06, -1.61], degrees=True).as_matrix()

With the rotation matrix ready, we need to calculate the offset and then we can perform the rotation.

In [11]:
# Calculate offset for rotation.
# out_center = rot_matrix @ ((np.array(data.shape) - 1) / 2)
# in_center = (np.array(data.shape) - 1) / 2
# offset = in_center - out_center

# Rotate data using spline interpolation as default.
# data = ndimage.affine_transform(data, rot_matrix, offset)

Note: Because scipy.ndimage.affine_transform uses spline interpolation and runs on a single CPU thread, rotating large volumes this way can take some time (hours).

Structure tensor

Now that we have the volume we would like to analyse using structure tensor in our data variable, there's a few important choices we need to make.

  1. Choose noise scale, $\sigma$ (sigma). Structures smaller than $\sigma$ will be removed by smoothing.
  2. Choose an integration scale, $\rho$ (rho). This is the size over the neighborhood in which the orientation is to be analysed.
  3. Choose data type/precision. 64-bit floating point precision is less likely to result in numerical issues, e.g. division by zero, but will use twice as much memory as 32-bit precision. As a rule of thumb, use 64-bit floats unless you're working with very large volumes and/or are running out of memory. To overcome numerical errors, one trick is to rescale the data or S before calculating eigenvectors.

Generally, float64 is favourable if you have sufficient memory and disk space, so we will use that here. Casting data from uint16 to float64 works fine for smaller volumes, however since it increases the since of the data four time, keeping everything in memory as we do here may cause memory issues for larger volumes.

In [12]:
data_f = data.astype(np.float64)

Chosing appropriate $\sigma$ and $\rho$ values is important for the results of the structure tensor analysis. In our case it makes sense to choose $\sigma$ based on the size of the fibers we want to analyze, so we set the fiber diameter and calculate $\sigma$ and $\rho$ based on that. For this type of data, choosing $\rho$ 2-4 times larger than $\sigma$ generally gives good results.

In [13]:
# Set fiber diameter.
fiber_diameter = 17 # μm

r = fiber_diameter / 2 / voxel_size
sigma = round(np.sqrt(r**2 / 2), 2)
rho = 4 * sigma

print('sigma:', sigma)
print('rho:', rho)

sigma: 0.74
rho: 2.96

It should be noted that chosing $\sigma$ much smaller than 0.5 may cause issues, as the filter kernel and its values, become very small. So we recommend keeping $\sigma > 0.5$, even if the structures are small.

Calculate S

Now the time has come to calculate the structure tensor, S. These calculations rely on the scipy.ndimage.gaussian_filter implementation and can be sped up significantly by moving the calculations to the GPU. The structure-tensor package supports this using CUDA through CuPy (see the end of this notebook).

In [14]:
%%time
S = structure_tensor_3d(data_f, sigma, rho)
print(f'Structure tensor information is carried in a {S.shape} array.')

Structure tensor information is carried in a (6, 500, 250, 596) array.
CPU times: user 23 s, sys: 3.02 s, total: 26 s
Wall time: 26.1 s

Having calculated S, the next step is to find the dominant structural orientations. To do so, we use the eig_special_3d function, which calculates the eigenvalues and eigenvectors of S. We set the argument full=False to get only the most important vector. Setting full=True will return a (3, 3, ...) array instead of a (3, ...) array instead.

In [15]:
%%time
val, vec = eig_special_3d(S, full=False)
print(f'Orientation information is carried in a {vec.shape} array.')

Orientation information is carried in a (3, 500, 250, 596) array.
CPU times: user 14.7 s, sys: 4.15 s, total: 18.8 s
Wall time: 18.9 s

Note that the first axis of the vector correspond to (z, y, x), meaning that vec[0] is the length of the vector along the Z-axis, not the X-axis as you might expect.

We will not be using valand S anymore here, so let's clear them to free up some memory.

In [16]:
# Remove variables to free up memory.
del val, S

Save eigenvectors

In case we wanted to view them externally, we could save the dominant vectors. Before we save the vectors we compress them to 8-bit precision to save disk space. If we want higher precision we can save them as 16-bit integer or floating point.

In [17]:
# temp_folder = '../tmp/'
# save_path = os.path.join(temp_folder, os.path.basename(file_path))

# if temp_folder and not os.path.exists(temp_folder):
#     os.mkdir(temp_folder)
In [18]:
# vec_uint8 = np.moveaxis(vec, 0, -1).copy()
# vec_uint8 *= 255
# vec_uint8 = vec_uint8.astype(np.uint8)

# ni_vol = nib.Nifti1Image(vec_uint8, affine=nii_file.affine)
# nib.save(ni_vol, save_path.replace('.nii', '-vec.nii'))

We could do something similar for the eigenvalues or other values we want to save to view externally or store for later use.

Orientation analysis

Now that we have the eigenvectors, we can use them to extract statistics about the orientation of the fibers, as well as to segment the different fiber classes. In our case we define fiber classes based on the direction of the bundles. The UD fibers, making up the majority, span along the X-axis (0 degree), while the three backing fiber classes are rotated -45, 45 and 90 degrees around the Z-axis, respectively.

We define four class vectors along the direction of the fibers they represent, in accordance with the format of vec, here (z, y, x).

In [19]:
# Names of the four fiber classes.
class_names = ['0', '45', '-45', '90']

# Unit vectors representing each of the four fiber classes.
class_vectors = np.array([[0, 0, 1], [0, 1, 1], [0, -1, 1], [0, 1, 0]], dtype=np.float64)
class_vectors /= np.linalg.norm(class_vectors, axis=-1)[..., np.newaxis]

for n, v in zip(class_names, class_vectors):
    print(f'Class vector ({n}) [z,y,x]: {v}')

Class vector (0) [z,y,x]: [0. 0. 1.]
Class vector (45) [z,y,x]: [0.         0.70710678 0.70710678]
Class vector (-45) [z,y,x]: [ 0.         -0.70710678  0.70710678]
Class vector (90) [z,y,x]: [0. 1. 0.]

With the class vectors defined, we call the calculate_angles helper function.

In [20]:
vec_class, eta_os, theta, out_xy_angle, in_xy_angle = calculate_angles(vec, class_vectors)

The function returns five metrics for each vector.

  1. vec_class contains a label for each voxel.
  2. eta_os contains a stiffness estimate, $\eta_O$, for each voxel
  3. theta contains the angle, $\theta$, between vec and class 0, at each voxel.
  4. out_xy_angle angle betweenvec and the XY-plane (out-of-plane orientation).
  5. in_xy_angle rotation of vec in the XY-plane (in-plane orientation).

Since we passed in the eigenvectors for the compelte volume, calculate_angles will give us the metrics for the entire volume. However, we don't have to do that. It is not necessary to calculate the metrics for all vectors at once. The function doesn't actually care about the shape of vec except for the first dimension, which must be three. This is because all the metrics are calculated individually for each vector. Thus, we can calculate the metrics for individual slices and visualize them.

First we grab some slices of data and their corresponding ST eigenvectors. Note: You can easily change/add/remove which slices to plot in the cell below.

In [21]:
data_slices = [
    data[data.shape[0] // 2].copy(),
    data[:, data.shape[1] // 2].copy(),
    data[:, :, data.shape[2] // 2].copy(),
    data[:, :, data.shape[2] // 3].copy(),
    data[:, :, data.shape[2] // 7].copy(),
]

vec_slices = [
    vec[:, vec.shape[1] // 2].astype(np.float32),
    vec[:, :, vec.shape[2] // 2].astype(np.float32),
    vec[:, :, :, vec.shape[3] // 2].astype(np.float32),
    vec[:, :, :, vec.shape[3] // 3].astype(np.float32),
    vec[:, :, :, vec.shape[3] // 7].astype(np.float32),
]

For each of the slices we calculate the metrics and draw them on top of the raw data. fig_with_colorbar is just a helper function for plotting, while show_metrics is responsible for calling calculate_angles for each slice and displaying it using fig_with_colorbar.

In [22]:
def fig_with_colorbar(d, o, title, alpha=0.5, cmap=None, vmin=None, vmax=None):
    """Creates a figure with data, overlay and color bar."""
    fig, ax = plt.subplots(1, 1, figsize=(10, 10))
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    ax.imshow(d, cmap='gray')
    if np.issubdtype(o.dtype, np.integer):
        cmap = plt.get_cmap('gist_rainbow', len(class_names))
        im = ax.imshow(o, alpha=alpha, cmap=cmap, vmin=-.5, vmax=len(class_names) - .5)
        cbar = fig.colorbar(im, cax=cax, orientation='vertical', ticks=np.arange(len(class_names)))
        cbar.ax.set_yticklabels(class_names)
    else:
        im = ax.imshow(o, alpha=alpha, cmap=cmap, vmin=vmin, vmax=vmax)
        fig.colorbar(im, cax=cax, orientation='vertical')
    ax.set_title(title)
    add_scalebar(ax, scale=voxel_size / 1000)
    plt.show()

def show_metrics(data_slices, vec_slices, fiber_threshold=None):
    
    for i, (d, v) in enumerate(zip(data_slices, vec_slices)):
        vec_class, eta_os, theta, out_xy_angle, in_xy_angle = calculate_angles(v, class_vectors)

        if fiber_threshold is not None:
            fiber_mask = d <= fiber_threshold

            vec_class = np.ma.masked_where(fiber_mask, vec_class)          
            eta_os[fiber_mask] = np.nan
            theta[fiber_mask] = np.nan
            out_xy_angle[fiber_mask] = np.nan
            in_xy_angle[fiber_mask] = np.nan

        fig_with_colorbar(d, vec_class, 'Class', alpha=0.7)
        fig_with_colorbar(d, eta_os, '$\eta_O$(X)', alpha=0.5, vmin=0, vmax=1)
        fig_with_colorbar(d, theta, 'Angle from X (0-90 deg.)', alpha=0.5, vmin=0, vmax=90)
        fig_with_colorbar(d, in_xy_angle, 'In-plane orientation (-90-90 deg.)', cmap='bwr', alpha=0.5, vmin=-90, vmax=90)
        fig_with_colorbar(d, out_xy_angle, 'Out-of-plane orientation (-90-90 deg.)', cmap='bwr', alpha=0.5, vmin=-90, vmax=90)
        fig_with_colorbar(d, theta, 'Angle from X (0-10 deg.)', alpha=0.5, vmin=0, vmax=10)
        fig_with_colorbar(d, in_xy_angle, 'In-plane orientation (-10-10 deg.)', cmap='bwr', alpha=0.5, vmin=-10, vmax=10)
        fig_with_colorbar(d, out_xy_angle, 'Out-of-plane orientation (-10-10 deg.)', cmap='bwr', alpha=0.5, vmin=-10, vmax=10)
        
show_metrics(data_slices, vec_slices)