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.
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:
../notebookscontains the notebooks.../originalsshould contain the original data files.../tmpshould contain any temperary files or output we create.
# 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)
Next step is to load the file using the nibabel package.
# 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')
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.
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.
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.
# 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.
# 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))
# Print the shape of the cut-out.
data_shape = data.shape
print('Cut-out shape:', data_shape)
Now that we have a smaller volume to work with, let's start by once again having a look at it.
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.
# 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.
# 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.
- Choose noise scale, $\sigma$ (
sigma). Structures smaller than $\sigma$ will be removed by smoothing. - Choose an integration scale, $\rho$ (
rho). This is the size over the neighborhood in which the orientation is to be analysed. - 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
Sbefore 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.
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.
# 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)
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).
%%time
S = structure_tensor_3d(data_f, sigma, rho)
print(f'Structure tensor information is carried in a {S.shape} array.')
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.
%%time
val, vec = eig_special_3d(S, full=False)
print(f'Orientation information is carried in a {vec.shape} array.')
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.
# 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.
# 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)
# 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).
# 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}')
With the class vectors defined, we call the calculate_angles helper function.
vec_class, eta_os, theta, out_xy_angle, in_xy_angle = calculate_angles(vec, class_vectors)
The function returns five metrics for each vector.
vec_classcontains a label for each voxel.eta_oscontains a stiffness estimate, $\eta_O$, for each voxelthetacontains the angle, $\theta$, betweenvecand class 0, at each voxel.out_xy_angleangle betweenvecand the XY-plane (out-of-plane orientation).in_xy_anglerotation ofvecin 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.
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.
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)
The figures above allow us that the method appears to work. Voxels are properly separated into the four classes and orientations seem correct. However, so far we've treated all voxels equally, but as we can see above, the volume is not just fibers. Classifying non-fiber voxels into one of the fiber classes or calculating their stiffness doesn't make sense and will only pollute our results.
Luckily, separating fibers from background can be done using a simple threshold.
Choosing fiber threshold¶
To choose a good intensity threshold for the fibers we can use a histogram and choose the value manually. Alternatively we can use a method such as Otsu's threshold. However, first let's have a look at the intensity distribution using a histogram.
# We have a lot of data, so let's have a reasonable large number of bins.
threshold_bins = 512
ax = plt.subplot(1, 1, 1)
ax.hist(data.flat, bins=threshold_bins, density=True)
plt.show()
Looking at the histogram there appears to be three different distributions. Let's have a look at the data again to see why this may be.
fig, ax = plt.subplots(1, 1, figsize=(10, 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()
From the image above we see that we have fiber, which is high intensity, non-fiber inside the sample, which is lower intensity and data outside the sample which has the lowest intensity. Since we want to separate the fibers from the rest (background), it would be sense to not include data from outside the sample, as the priority is to separate the two distributions inside the sample.
Let's crop away the edges of our data see if this removes the low intensity distribution. We will put this subset of data in the threshold_data variable, which is only used for choosing the threshold.
# Crop away data close to the edges that may be outside the fiber material.
# The two distributions we want to separate are fiber and non-fiber "inside" the material.
threshold_data = data[:, :, 35:-25]
fig, (ax_h, ax_im) = plt.subplots(1, 2, figsize=(10, 3))
ax_h.hist(threshold_data.flat, bins=threshold_bins, density=True)
i = threshold_data.shape[0] // 2
img = threshold_data[i]
ax_im.imshow(img.squeeze(), cmap='gray')
ax_im.set_title('YZ-slice')
ax_im.set_xlabel('Z')
ax_im.set_ylabel('Y')
add_scalebar(ax_im, scale=voxel_size / 1000, ticks=True)
plt.show()
Now that we've selected some data that is strictly inside the sample, the low intensity distribution is gone and we just have to find the best separation between what looks like two distributions. For this we can use Otsu's threshold and choose one ourselves.
otsu_threshold = filters.threshold_otsu(threshold_data.reshape(-1, 1), nbins=threshold_bins)
hand_picked_threshold = 41000
print('Otsu threshold:', otsu_threshold)
print('Hand-picked threshold:', hand_picked_threshold)
The we plot a histogram for the threshold data with both thresholds marked.
ax = plt.subplot(1, 1, 1)
ax.hist(threshold_data.flat, bins=threshold_bins, density=True)
ax.axvline(otsu_threshold, c='r', label='Otsu')
ax.axvline(hand_picked_threshold, c='g', label='Hand-picked')
plt.legend()
plt.show()
Fiber fraction¶
We can also easily calculate the fiber fraction of the material, which may help us pick an appropriate threshold.
print('Otsu fiber fraction:', round(np.count_nonzero(threshold_data >= otsu_threshold) / threshold_data.size, 3))
print('Hand-picked fiber fraction:', round(np.count_nonzero(threshold_data >= hand_picked_threshold) / threshold_data.size, 3))
Here we will be using our hand-picked threshold. Let's have a look at the threshold mask.
fiber_threshold = hand_picked_threshold
threshold_mask = data >= fiber_threshold
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
i = data_shape[2] // 2
ax.imshow(data[i], cmap='gray')
ax.imshow(threshold_mask[i], alpha=0.2, cmap='cool')
ax.set_title('YZ-slice')
ax.set_xlabel('Z')
ax.set_ylabel('Y')
add_scalebar(ax, length=1, scale=voxel_size / 1000, ticks=True)
plt.show()
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
i = data_shape[1] // 2
ax.imshow(data[:, i], cmap='gray')
ax.imshow(threshold_mask[:, i], alpha=0.2, cmap='cool')
ax.set_title('XZ-slice')
ax.set_xlabel('Z')
ax.set_ylabel('X')
add_scalebar(ax, length=1, scale=voxel_size / 1000, ticks=True)
plt.show()
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
i = data_shape[0] // 2
ax.imshow(data[..., i], cmap='gray')
ax.imshow(threshold_mask[..., i], alpha=0.2, cmap='cool')
ax.set_title('XY-slice')
ax.set_xlabel('Y')
ax.set_ylabel('X')
add_scalebar(ax, length=1, scale=voxel_size / 1000, ticks=True)
plt.show()
Now, we can show the figures from before, but without coloring background, which gives a much clearer view of the fiber structures and their orientation. In practice we have created another class, background, in addition to the four fiber classes.
show_metrics(data_slices, vec_slices, fiber_threshold=hand_picked_threshold)
Histograms¶
The figures above are great for validating that we've done the calculations correctly and that our results are meaningful. However, even with this many figures we're only actually showing a very small fraction of the volume. To get a proper overview of the fiber orientations and classes we need to aggregate the information in a meaningful way. One way to do this is to use histograms.
By combining the threshold mask (indicating what is fiber) with the class labels, we can create a fiber class mask for each class.
# Count number of fiber voxels.
fiber_count = np.count_nonzero(threshold_mask)
# Create a mask for each fiber class.
vec_class_masks = [threshold_mask] + [(vec_class == l) & threshold_mask for l in range(len(class_vectors))]
Fiber fractions¶
Using these masks we calculate the fraction of the fibers belonging to each class.
fiber_fractions = [np.count_nonzero(m) / fiber_count for m in vec_class_masks[1:]]
for i, f in enumerate(fiber_fractions):
c = class_names[i]
print(f'Fraction ({c}): {round(f, 6)}')
plt.bar(class_names, fiber_fractions, align='center', width=0.5)
plt.xticks(range(len(fiber_fractions)))
plt.show()
Angle from X ($\theta$)¶
The histograms below show the distribution of the absolute angles between the fibers and the X-axis unit vector, along with the median angle. The resoulution of the histograms is determined by bins_per_degree.
# Four bins per degree = 0.25 degree resolution.
bins_per_degree = 4
# Calculate theta for each class using the masks.
thetas = [theta[m] for m in vec_class_masks]
def hist(i, ax):
if i == 0:
title = 'All classes'
ax.set_xlim([0, 90])
else:
class_name = class_names[i - 1]
title = f'Class {class_name}'
class_value = int(class_name)
ax.set_xlim([max(0, abs(class_value) - 30) , min(90, abs(class_value) + 30)])
m = np.median(thetas[i])
ax.axvline(m, c='b', ls='-', label=f'$Med={round(m, 2)}\degree$')
plt.legend()
ax.set_title(title)
ax.hist(thetas[i], density=True, bins=np.arange(0, 90 + 1 / bins_per_degree, 1 / bins_per_degree))
ax.set_xlabel('Angle (deg.)')
ax.set_ylabel('Fraction')
fig = plt.figure(figsize=(15, 7))
fig.suptitle('Angle from X')
for i in range(len(class_names) + 1):
ax = plt.subplot(2, 3, i + 1)
hist(i, ax)
plt.tight_layout()
plt.subplots_adjust(top=0.85)
plt.show()
In-plane orientation (rotaiton in XY-plane)¶
The histograms below show the distribution of the rotation around the Z-axis (in-XY-plane rotation), the mean and the standard deviation for each class.
in_xy_angles = [in_xy_angle[m] for m in vec_class_masks]
def hist(i, ax):
if i == 0:
title = 'All classes'
class_value = 0
ax.set_xlim([-90, 90])
ax.axvline(0, c='k', ls='--')
ax.axvline(-45, c='k', ls='--')
ax.axvline(45, c='k', ls='--')
else:
class_name = class_names[i - 1]
title = f'Class {class_name}'
class_value = int(class_name)
ax.set_xlim([class_value - 20, class_value + 20])
ax.axvline(class_value, c='k', ls='--', label=f'${class_value}\degree$')
if class_value == 90:
ax.set_xticks(range(class_value - 30, class_value + 31, 10))
ax.set_xticklabels([60, 70, 80, '$\pm90$', -80, -70, -60])
bins = np.arange(0, 180 + 1 / bins_per_degree, 1 / bins_per_degree)
values = in_xy_angles[i].copy()
values[values < 0] = 90 + (values[values < 0] + 90)
else:
bins = np.arange(-90, 90 + 1 / bins_per_degree, 1 / bins_per_degree)
values = in_xy_angles[i]
ax.set_title(title)
ax.hist(values, density=True, bins=bins)
mean = np.mean(values)
std = np.std(values)
ax.set_xlabel('Angle (deg.)')
ax.set_ylabel('Fraction')
if i != 0:
ax.axvline(mean, c='r', ls='-', label=f'$\\bar{{x}}={round(mean if mean < 90 else -180 + mean, 2)}\degree$')
ax.axvline(mean - std, c='r', ls='--', label=f'$s=\pm{round(std, 2)}\degree$')
ax.axvline(mean + std, c='r', ls='--')
plt.legend()
fig = plt.figure(figsize=(15, 7))
fig.suptitle('In-plane orientation')
for i in range(len(class_names) + 1):
ax = plt.subplot(2, 3, i + 1)
hist(i, ax)
plt.tight_layout()
plt.subplots_adjust(top=0.85)
plt.show()
Out-of-plane orientation (angle from xy-plane)¶
The histograms below show the distribution of the angle from the out-of-plane (XY-plane) angle.
out_xy_angles = [out_xy_angle[m] for m in vec_class_masks]
def hist(i, ax):
if i == 0:
title = 'All classes'
class_value = None
else:
title = f'Class {class_names[i - 1]}'
class_value = int(class_names[i - 1])
ax.set_title(title)
mean = np.mean(out_xy_angles[i])
std = np.std(out_xy_angles[i])
if class_value == 90:
ax.set_xlim([-20, 20])
else:
ax.set_xlim([-10, 10])
ax.hist(out_xy_angles[i], density=True, bins=np.arange(-90, 90 + 1 / bins_per_degree, 1 / bins_per_degree))
ax.set_xlabel('Angle (deg.)')
ax.set_ylabel('Fraction')
ax.axvline(0, c='k', ls='--', label='$0\degree$')
if i != 0:
ax.axvline(mean, c='r', ls='-', label=f'$\\bar{{x}}={round(mean, 2)}\degree$')
ax.axvline(mean - std, c='r', ls='--', label=f'$s=\pm{round(std, 2)}\degree$')
ax.axvline(mean + std, c='r', ls='--')
plt.legend()
fig = plt.figure(figsize=(15, 7))
fig.suptitle('Out-of-plane orientation')
for i in range(len(class_names) + 1):
ax = plt.subplot(2, 3, i + 1)
hist(i, ax)
plt.tight_layout()
plt.subplots_adjust(top=0.85)
plt.show()
The class fractions and orientation distributions are useful for determining if the material follows production specifications. Here we've chosen a specific set of metrics to calculate and plot, but many other distributions can just as easily be obtained and plotted.
Handling larger volumes¶
So far we've been using a small volume to demonstrate. As we've mentioned, the reason is that the calculations and their results are quite memory intensive. One way to get around this is to separate the volume into blocks and handle them separately.
For this approach to work, we need to pad each block (in the code called "crop") depending on the filter kernel size, which is determined by sigma and rho. To help with this the structure_tensor_workers module implements a number of functions. The code below demonstrates how these functions can be used to generate blocks, calculate the eigenvectors for the block and insert them into a complete vector volume.
from structure_tensor_workers import get_crop_generator, insert_crop, remove_padding
We create a new vec array large enough to contain 3D vectors for the entire volume. We use uint8 to limit the memory usage of the array. This should be fine if we're just using the vectors for visuals.
class_labels = np.zeros(nii_file.shape, dtype=np.uint8)
print('Shape:', class_labels.shape)
print('Bytes:', class_labels.nbytes)
We new use the get_crop_generator function to create a generator, which will generate all the blocks for the nii_file.dataobj. The block size, crop_size, determines the size of the block, excluding the padding. Larger blocks will result in less overhead from padding, but higher memory usage during calculations.
We can plot the blocks to better get an understanding of how they are created.
crop_size = 150
def plot_cube(i, mode, ax, bbox, surface=False):
bbox = np.asarray(bbox)
mins = bbox.min(axis=0)
maxs = bbox.max(axis=0)
points = np.array(np.meshgrid([mins[0], maxs[0]], [mins[1], maxs[1]], [mins[2], maxs[2]])).T.reshape(-1,3)
edges = [
[points[0], points[1], points[5], points[4]],
[points[0], points[2], points[6], points[4]],
[points[0], points[1], points[3], points[2]],
[points[2], points[3], points[7], points[6]],
[points[1], points[3], points[7], points[5]],
[points[4], points[5], points[7], points[6]],
]
c = plt.cm.Set1(i)[:-1]
if mode == 'vol':
faces = Poly3DCollection(edges, linewidths=3, edgecolor='k')
faces.set_facecolor((1, 1, 1, 0.01))
if surface:
# Plot back surface.
xx, zz = np.meshgrid(np.arange(mins[0], maxs[0]), np.arange(mins[2], maxs[2]))
X = xx
Z = zz
Y = np.full(X.shape, maxs[1])
d = data[:, maxs[1] - 1].T - np.min(data)
d = d / np.max(d)
color = plt.cm.gray(d)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=color, shade=False)
# Plot YZ-surface.
yy, zz = np.meshgrid(np.arange(mins[1], maxs[1]), np.arange(mins[2], maxs[2]))
Y = yy
Z = zz
X = np.full(Y.shape, (mins[0] + maxs[0]) // 2)
d = data[X.flat[0]].T - np.min(data)
d = d / np.max(d)
color = plt.cm.gray(d)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=color, shade=False)
# Plot XY-surface.
# xx, yy = np.meshgrid(np.arange(mins[0], maxs[0]), np.arange(mins[1], maxs[1]))
# X = xx
# Y = yy
# Z = np.full(X.shape, (mins[2] + maxs[2]) // 2)
# d = data[..., Z.flat[0]].T - np.min(data)
# d = d / np.max(d)
# color = plt.cm.gray(d)
# ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=color, shade=False)
elif mode == 'pad':
faces = Poly3DCollection(edges, linewidths=1, ls=':', edgecolor=c + (1 / (i + 1),))
faces.set_facecolor(c + (0.1 / (i + 1),))
else:
faces = Poly3DCollection(edges, linewidths=2, edgecolor=c + (1 / (i + 1),))
faces.set_facecolor(c + (0.2 / (i + 1),))
ax.add_collection3d(faces)
# Plot the points themselves to force the scaling of the axes
ax.scatter(points[:,0], points[:,1], points[:,2], s=0)
fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(1, 1, 1, projection='3d', proj_type='ortho')
ax.view_init(elev=10, azim=-75)
plt.axis('off')
# Plot volume.
# draw_vol = nii_file.dataobj[:128, :128, :128]
draw_vol = data
plot_cube(-1, 'vol', ax, [(0,0,0), draw_vol.shape], surface=True)
# Plot some blocks.
for i, (crop, pos, pad) in enumerate(get_crop_generator(draw_vol, max(sigma, rho), crop_size=crop_size)):
if i > 0:
break
plot_cube(i, 'pos', ax, np.minimum(np.array(pos).T, np.array(draw_vol.shape)[np.newaxis]))
pad[:, 0] = -pad[:, 0]
plot_cube(i, 'pad', ax, np.minimum(np.array(pos + pad).T, np.array(draw_vol.shape)[np.newaxis]))
ax.set_ylim(ax.get_xlim())
ax.set_zlim(ax.get_xlim())
plt.savefig('figures/blocks1_surf.png', bbox_inches='tight', pad_inches=0)
plt.show()
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, projection='3d', proj_type='ortho')
ax.view_init(elev=10, azim=-75)
plt.axis('off')
# Plot volume.
# draw_vol = nii_file.dataobj[:128, :128, :128]
draw_vol = data
plot_cube(-1, 'vol', ax, [(0,0,0), draw_vol.shape])
# Plot some blocks.
for i, (crop, pos, pad) in enumerate(get_crop_generator(draw_vol, max(sigma, rho), crop_size=crop_size)):
if i > 0:
break
plot_cube(i, 'pos', ax, np.minimum(np.array(pos).T, np.array(draw_vol.shape)[np.newaxis]))
pad[:, 0] = -pad[:, 0]
plot_cube(i, 'pad', ax, np.minimum(np.array(pos + pad).T, np.array(draw_vol.shape)[np.newaxis]))
ax.set_xlim(np.asarray(ax.get_xlim()) / 4)
ax.set_ylim(ax.get_xlim())
ax.set_zlim(ax.get_xlim())
plt.savefig('figures/blocks1.pdf', bbox_inches='tight', pad_inches=0)
plt.show()
The inner red cube is a volume of the size specified by crop_size. The values inside the inner red cube are the ones that will be valid and that we can use afterward. However, the crop volume returned by get_crop_generator is actually larger and covers the outer red box. The difference between them is calculated from max(sigma, rho) og the truncate argument of the scipy.ndimage.gaussian_filter function, which is 4 by default. This is because the kernel_radius used in the filter is calculated as kernel_radius = int(sigma * truncate + 0.5).
truncate = 4
kernel_radius = int(max(sigma, rho) * truncate + 0.5)
print('Padding needed:', kernel_radius)
To calculate S for a voxel the filters use 12 neighboring voxels in each of the six directions. These neighbors should be the same no matter the block size. To ensure this we include the extra voxels when calculating S, however the values for the padding voxels may not be correct and we have to throw them away. As can be seen in the figure we don't pad on the edge of the volume where there's actually no data. Padding here is handled by the scipy.ndimage.gaussian_filter function. We only pad the blocks with actual data.
Let's draw some more blocks.
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, projection='3d', proj_type='ortho')
ax.view_init(elev=10, azim=-75)
plt.axis('off')
# Plot volume.
# draw_vol = nii_file.dataobj[:128, :128, :128]
draw_vol = data
plot_cube(-1, 'vol', ax, [(0,0,0), draw_vol.shape])
# Plot some blocks.
for i, (crop, pos, pad) in enumerate(get_crop_generator(draw_vol, max(sigma, rho), crop_size=crop_size)):
if i > 1:
break
plot_cube(i, 'pos', ax, np.minimum(np.array(pos).T, np.array(draw_vol.shape)[np.newaxis]))
pad[:, 0] = -pad[:, 0]
plot_cube(i, 'pad', ax, np.minimum(np.array(pos + pad).T, np.array(draw_vol.shape)[np.newaxis]))
ax.set_xlim(np.asarray(ax.get_xlim()) / 2)
ax.set_ylim(ax.get_xlim())
ax.set_zlim(ax.get_xlim())
plt.savefig('figures/blocks2.pdf', bbox_inches='tight', pad_inches=0)
plt.show()
... and some more.
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, projection='3d', proj_type='ortho')
ax.view_init(elev=10, azim=-75)
plt.axis('off')
# Plot volume.
# draw_vol = nii_file.dataobj[:512, :512, :512]
draw_vol = data
plot_cube(-1, 'vol', ax, [(0,0,0), draw_vol.shape])
# Plot some blocks.
for i, (crop, pos, pad) in enumerate(get_crop_generator(draw_vol, max(sigma, rho), crop_size=crop_size)):
if i > 30:
break
plot_cube(i, 'pos', ax, np.minimum(np.array(pos).T, np.array(draw_vol.shape)[np.newaxis]))
pad[:, 0] = -pad[:, 0]
plot_cube(i, 'pad', ax, np.minimum(np.array(pos + pad).T, np.array(draw_vol.shape)[np.newaxis]))
ax.set_ylim(ax.get_xlim())
ax.set_zlim(ax.get_xlim())
plt.savefig('figures/blocks3.pdf', bbox_inches='tight', pad_inches=0)
plt.show()
Once the eigenvectors have been caculated for the block, we rescale them and change their type to uint8 to save space. The we use the fiber_threshold to create a mask, which we will use to only insert vector values into vec where there's fiber material. After that we remove the padding and insert the eigenvectors into vec at the correct location using the insert_crop function.
Here we will only run the 30 first blocks, as doing them all takes a long time.
%%time
count = 0
for crop, pos, pad in get_crop_generator(nii_file.dataobj, max(sigma, rho), crop_size=crop_size):
# Rorate crop.
crop = np.rot90(crop, k=3, axes=(0, 2))
# Calculate structure tensor.
S_crop = structure_tensor_3d(crop.astype(np.float64), sigma, rho)
val_crop, vec_crop = eig_special_3d(S_crop)
# Free memory.
del val_crop, S_crop
# Calculate metrics.
vec_class, eta_os, theta, out_xy_angle, in_xy_angle = calculate_angles(vec_crop, class_vectors)
# Create mask.
mask = crop >= fiber_threshold
# Set background.
vec_class += 1
vec_class[mask] = 0
# Rotate back.
vec_class = np.rot90(vec_class, k=1, axes=(0, 2))
# Insert class labels.
insert_crop(class_labels, vec_class, pos, pad)
count += 1
if count >= 30:
break
print(f'Volume with size {class_labels.shape} was split into {count} blocks.')
Because we stop early, only a small part of the volume has been segmented. To run it all, simply remove the break statement.
cl = class_labels[..., 150].copy()
cl -= 1
cl = np.ma.masked_where(cl == 255, cl)
fig_with_colorbar(nii_file.dataobj[..., 150], cl, 'Class', alpha=0.3)
CUDA support¶
One way to speed up computations is by using GPUs. The structure-tensor package support using CUDA through CuPy. This requires an Nvidia GPU and CUDA to be installed on your system. If you have a powerfull GPU this should speed up calculations significantly, but you have to be carefull with memory, as GPUs generally has less memory than your CPU.
Below we demonstrate how the same calculations could be done using functions from the structure_tensor.cp module.
try:
import cupy as cp
from structure_tensor.cp import eig_special_3d as eig_special_3d_cp, structure_tensor_3d as structure_tensor_3d_cp
has_cupy = True
except:
print('CuPy import failed')
has_cupy = False
if has_cupy:
class_labels = np.zeros(nii_file.shape, dtype=np.uint8)
print('Shape:', class_labels.shape)
print('Bytes:', class_labels.nbytes)
%%time
if has_cupy:
count = 0
for crop, pos, pad in get_crop_generator(nii_file.dataobj, max(sigma, rho), crop_size=crop_size):
# Rorate crop.
crop = np.rot90(crop, k=3, axes=(0, 2))
# Calculate structure tensor.
S_crop = structure_tensor_3d_cp(crop.astype(np.float64), sigma, rho)
val_crop, vec_crop = eig_special_3d_cp(S_crop)
# Free memory.
del val_crop, S_crop
# Move to NumPy.
vec_crop = cp.asnumpy(vec_crop)
# Calculate metrics.
vec_class, eta_os, theta, out_xy_angle, in_xy_angle = calculate_angles(vec_crop, class_vectors)
# Create mask.
mask = crop >= fiber_threshold
# Set background.
vec_class += 1
vec_class[mask] = 0
# Rotate back.
vec_class = np.rot90(vec_class, k=1, axes=(0, 2))
# Insert class labels.
insert_crop(class_labels, vec_class, pos, pad)
count += 1
if count >= 30:
break
print(f'Volume with size {class_labels.shape} was split into {count} blocks.')
else:
print('CuPy not availble.')
cl = class_labels[..., 150].copy()
cl -= 1
cl = np.ma.masked_where(cl == 255, cl)
fig_with_colorbar(nii_file.dataobj[..., 150], cl, 'Class', alpha=0.3)
By combining blocking with parallism and CUDA, we can speed up calculations, so even large volumes can be analysed in a few minutes. Our other notebook demonstrates how this can be done.