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.
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.
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
.
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
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.# 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.
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()
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()
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).
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.
sigma
). Structures smaller than $\sigma$ will be removed by smoothing.rho
). This is the size over the neighborhood in which the orientation is to be analysed.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.
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.
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 val
and S
anymore here, so let's clear them to free up some memory.
# Remove variables to free up memory.
del val, S
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.
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_class
contains a label for each voxel. eta_os
contains a stiffness estimate, $\eta_O$, for each voxeltheta
contains the angle, $\theta$, between vec
and class 0, at each voxel.out_xy_angle
angle betweenvec
and the XY-plane (out-of-plane orientation).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.
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.
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()
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)
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))]
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()
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()
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()
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.
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)
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.