Analysis of HF401TT-13_FoV16.5_Stitch.txm¶
Author: Niels Jeppesen (niejep@dtu.dk)
This notebook can be used to reproduce the results in the Characterization of the Fiber Orientations in Non-Crimp Glass Fiber Reinforced Composites using Structure Tensor paper. The notebook is based on the StructureTensorFiberAnalysisAdvancedDemo notebook, please see that for more details.
This notebook reads the TXM data format. In our experience NIfTI is much faster to read than TXM using the packages we use here.
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 tqdm structure-tensor. To start the notebook we will of course also need jupyter. Install the dxchange package from GitHub, see note below.
Note: The current release of the dxchange package has a bug in the TXM reader. This is fixed in the GitHub master branch, so for now the package must be installed directly from GitHub using pip install git+https://github.com/data-exchange/dxchange.git. The setup.py does not specify requirements for the package, which can be found in dxchange requirements.txt instead.
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
from multiprocessing import Pool
import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import numpy as np
from dxchange import reader
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from scipy import ndimage
from scipy.spatial.transform import Rotation as R
from skimage import filters
from structure_tensor import eig_special_3d, structure_tensor_3d
from tqdm import tqdm
from structure_tensor_workers import calculate_angles, get_crops, init_worker, structure_tensor_analysis_v1
Load and display data¶
First, we'll load a sample of the data and some meta data. 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.../notebooks/figures/<file_name>contains optionally saved figures.
# Set file name and path.
file_name = 'HF401TT-13_FoV16.5_Stitch.txm'
file_path = os.path.join('../originals/', file_name)
# Change this to True to save figures to folder.
save_figures = False
# Create folder for to save figures.
fig_path = os.path.join('figures', os.path.basename(file_name))
if save_figures and fig_path and not os.path.exists(fig_path):
os.makedirs(fig_path)
# Read the first slice and metadata.
sample, meta = reader.read_txm(file_path, slice_range=(1, None, None))
data_shape = (meta['number_of_images'], meta['image_height'], meta['image_width'])
data_type = sample.dtype
voxel_size = meta['pixel_size']
fiber_diameter = 17
print('Shape:', data_shape)
print('Data type:', data_type)
print('Voxel size:', voxel_size, 'μm')
print('Fiber diameter size:', fiber_diameter, 'μm')
Let's have a need look at the data, to ensure it was loaded correctly.
fig, axs = plt.subplots(1, 2, figsize=(20, 20))
i = data_shape[0] // 3
img, _ = reader.read_txm(file_path, slice_range=((i, i + 1), None, None))
axs[0].imshow(img.squeeze(), cmap='gray')
i = data_shape[0] // 3 * 2
img, _ = reader.read_txm(file_path, slice_range=((i, i + 1), None, None))
axs[1].imshow(img.squeeze(), cmap='gray')
plt.show()
Crop data¶
Clearly a large part of the data is not of interest to us, so let's crop out the area of interest. You can cut out any part of the data you want to analyze.
x_slice = slice(200, data_shape[0] - 200)
y_slice = slice(80, data_shape[1] - 80)
z_slice = slice(700, data_shape[2] - 700)
To cut out and examine the two outermost backing layer individually, use the uncomment the slice setting below.
# Backing 1 (upper)
# x_slice = slice(200, data_shape[0] - 200)
# y_slice = slice(80, data_shape[1] - 80)
# z_slice = slice(data_shape[2] - 850, data_shape[2] - 700)
# Backing 2 (lower)
# x_slice = slice(200, data_shape[0] - 200)
# y_slice = slice(80, data_shape[1] - 80)
# z_slice = slice(700, 850)
Let's have a look at the cropped area. The TXM reader can be a bit slow, so we will only look at YZ-slices for now. We may have to come back and adjust the cropping on the X-axis later.
fig, axs = plt.subplots(1, 2, figsize=(20, 20))
i = data_shape[0] // 3
img, _ = reader.read_txm(file_path, slice_range=((i, i + 1), y_slice, z_slice))
axs[0].imshow(img.squeeze(), cmap='gray')
i = data_shape[0] // 3 * 2
img, _ = reader.read_txm(file_path, slice_range=((i, i + 1), y_slice, z_slice))
axs[1].imshow(img.squeeze(), cmap='gray')
plt.show()
Create new shape. Note: We rotate 180 degrees around the X axes to match convention.
data_shape = (x_slice.stop - x_slice.start, y_slice.stop - y_slice.start, z_slice.stop - z_slice.start)
print('New shape:', data_shape)
Save as RAW¶
To work with the data, we will save it as a raw data file, which will allow us to access the data using numpy.memmap. Using a memory-mapped file is both convenient and fast, especially since we will be accessing the data many times from different processes.
We will be placing all our processed data in a temporary folder, as it can be re-created using this notebook if deleted. We will also use data in the temp folder as a cache, so we don't have to create the RAW volume if it already exists. Thus, if you change the volume cropping etc., be sure to delete the RAW volume from the temp folder and create a new one.
temp_folder = '../tmp/'
if temp_folder and not os.path.exists(temp_folder):
os.mkdir(temp_folder)
data_path = os.path.join(temp_folder, file_name + f'_{x_slice.start}-{x_slice.stop}_{y_slice.start}-{y_slice.stop}_{z_slice.start}-{z_slice.stop}.raw')
if not os.path.exists(data_path):
# Number of images read at a time.
# Reduce this to lower memory usage.
chunck_size = 1024
print('Creating new file:', data_path)
data = np.memmap(data_path, dtype=data_type, shape=data_shape, mode='w+')
for i in tqdm(range(0, data_shape[0], chunck_size)):
x0 = i + x_slice.start
x1 = min(x_slice.stop, x0 + chunck_size)
chunck, _ = reader.read_txm(file_path, slice_range=((x0, x1), y_slice, z_slice))
data[i:i + chunck_size] = np.rot90(chunck, k=2, axes=(1, 2))
else:
print('Using existing file:', data_path)
data = np.memmap(data_path, dtype=data_type, shape=data_shape, mode='r')
Let's have a look some slices from all three axes.
Important: Remember to look at figures below to verify that the correct data is loaded into the data variable. If you've changed the slicing or the data you're working with, you may want to delete the RAW data cached in the temp_folder folder.
plt.figure(figsize=(20, 20))
ax = plt.subplot(2, 2, 1)
ax.imshow(data[data.shape[0] // 3], cmap='gray')
ax = plt.subplot(2, 2, 2)
ax.imshow(data[data.shape[0] // 3 * 2], cmap='gray')
ax = plt.subplot(2, 2, 3)
ax.imshow(data[:, data.shape[1] // 2], cmap='gray')
ax = plt.subplot(2, 2, 4)
ax.imshow(data[..., data.shape[-1] // 2], cmap='gray')
plt.show()
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. For choosing the threshold we make sure only to include inside-sample data. This makes it easier to separate foreground (fibers) and background inside the sample.
threshold_bins = 512
threshold_data = data[:, 50:-50, 50:-50]
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)
ax = plt.subplot(1, 1, 1)
ax.hist(threshold_data.flat, bins=threshold_bins, density=True)
ax.axvline(otsu_threshold, c='r')
ax.axvline(hand_picked_threshold, c='g')
plt.show()
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))
Structure tensor¶
We set $\sigma$ and $\rho$ based on the size of the fibers that we want to analyze. For more details see the related paper and StructureTensorFiberAnalysisDemo notebook.
Then, we set the block size (crop_size), which will determine the maximum size of the blocks we will partition the data into for the ST calculation. The maximum black size will depend on crop_size, the truncate value used for Gaussian filters (default is 4), and $\rho$ or $\sigma$ (usually $\rho$ because ti is largest).
We also set the fiber_threshold and specify a list of devices to use for the calculation. The list determines how the blocks will be distributed between devices for calculations. If we have a dual GPU system supporting CUDA, we can specify ['cuda:0', 'cuda:1'], which will distribute blocks evenly between the two GPUs. ['cuda:0', 'cuda:1', 'cuda:1'] would move two thirds to GPU 1, while ['cuda:0', 'cuda:1', 'cpu'] will move one third of the blocks to the CPU. Remember to update the device list to match your system resources. Specifying CUDA devices that are not available will result in exceptions and/or undefined behaviour. If you don't have CUDA GPUs available, just set device = ['cpu'] and specify the number of processes later.
The class_vectors area used to segment the voxels based on the orientations of the ST eigenvectors. It is a list of unit vectors, which represent each of the fiber classes. The first vector is considered the primary orientation and will be used for calculating orientation metrics for each voxel. For more details see the related paper and StructureTensorFiberAnalysisDemo notebook.
Lastly, bins_per_degree determines the number of bins used per degree when aggregating the orientation metrics.
r = fiber_diameter / 2 / voxel_size
sigma = round(np.sqrt(r**2 / 2), 2)
rho = 4 * sigma
print('sigma:', sigma)
print('rho:', rho)
crop_size = 230
print('crop_size:', crop_size)
print('Maximum block size:', crop_size + int(4 * max(rho, sigma) + 0.5))
fiber_threshold = hand_picked_threshold
# Important: Listing invalid CUDA devices may results in exceptions.
device = ['cuda:0', 'cuda:1', 'cuda:2', 'cuda:3']
class_names = ['0', '45', '-45', '90']
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]
bins_per_degree = 4
We can use the get_crops function to partition the data into blocks (crops) with proper padding. However, we will actually only use len(crops) here, as we will use multiprocessing to distribute the blocks accross multiple devices. We may include a function just for calculating the number of blocks/crops at a later point.
# Get crops as memory views.
crops, crop_positions, crop_paddings = get_crops(data, max(sigma, rho), crop_size=crop_size)
We will be using the structure_tensor_analysis_v1 function to calculate the structure tensor, $S$, and do the eigendecomposition for the blocks in parallel. Furthermore the function will aggregate the metrics:
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).
For this it uses the calculate_angles function, which is explained further in the StructureTensorFiberAnalysisDemo notebook. However, the metrics returned by structure_tensor_analysis_v1 have been aggregated (binned), as returning all the metrics for each block unaggregated would obvously consume large amounts of memory and be infeasible for large volumes. We will combine the aggregated data from each block afterwards.
However, it is possible to save the actual metrics for each voxel to the disk using memory maps as shown below. In the code below we create memory-mapped files for the eigenvectors, along with three other metrics. The structure_tensor_analysis_v1 function will be using these to save the metrics for the volume straight to the disk, which may require a significant amount of disk space, but shouldn't result in memory issues. See structure_tensor_workers.py for details. Saving the metrics to disk is optional. If you don't need the per voxel orientations for later, don't create the memory-maps and remove the entries from the init_args dictionary below. This will save you a lot of disk space and probably reduce the processing time.
# Output names, dtypes and shapes.
map_names = ['vec']
map_dtypes = [np.float16]
map_shapes = [(3, ) + data.shape]
# Output paths.
map_paths = {n: data_path.replace('.raw', f'-{n}-{sigma}-{rho}-{fiber_threshold}.raw') for n in map_names}
map_paths['vec'] = data_path.replace('.raw', f'-vec-{sigma}-{rho}.raw')
# Create maps.
maps = {}
new_maps = {}
for n, dtype, shape in zip(map_names, map_dtypes, map_shapes):
path = map_paths[n]
create_map = not os.path.exists(path) or os.stat(path).st_size != np.product(shape) * np.dtype(dtype).itemsize
if create_map:
print(f'Creating memory mapped file:\n{path}')
mmap = np.memmap(path, dtype=dtype, shape=shape, mode='w+')
if np.issubdtype(dtype, np.floating) and n not in ['S', 'val', 'vec']:
mmap[:] = np.nan
new_maps[n] = (path, dtype)
maps[n] = np.memmap(path, dtype=dtype, shape=shape, mode='r')
# Get vec memory mapped array.
vec = maps['vec']
Now we're finally ready to perform the analysis. We will be using multiprocessing.Pool to distribute the work across multiple CPU cores. If specified in the device list, the CPU will offload the majority of the work to a GPU, otherwise it'll do the calculations itself. Here, we will create four processes for each device in the list (processes=4 * len(device)). As our device list contains four GPUs, we will be starting 16 processes, four for each GPU. Beware that this may require a significant amount of GPU memory. If you experience out of memory exceptions, either reduce crop_size and/or the number of processes per device. Change the number of processes to fit your system.
if __name__ == '__main__':
results = []
init_args = {'data': data_path,
'dtype': data.dtype,
'shape': data.shape,
'rho': rho,
'sigma': sigma,
'crop_size': crop_size,
'threshold': fiber_threshold,
'class_vectors': class_vectors,
'bins_per_degree': bins_per_degree,
'device': device,
'return_aggregate': True
}
# Add output maps to args.
for n in new_maps:
path, dtype = new_maps[n]
init_args[n] = path
init_args[f'{n}_dtype'] = dtype
with Pool(processes=4 * len(device), initializer=init_worker, initargs=(init_args,)) as pool:
for res in tqdm(pool.imap_unordered(structure_tensor_analysis_v1, range(len(crops)), chunksize=1), total=len(crops)):
results.append(res)
All shared variables are passed to the workers using the init_args dictionary and init_workerfunction. This function is only called once per worker. When called, it will create all the memory-maps needed to read and write data during calculations, based on init_args. We will use the pool.imap_unordered function to tell the workers to perform the calculations for a block, by passing the structure_tensor_analysis_v1 function along with the index of the block which we want to analyze. The returned results are a tuple of aggregated metrics, as decribed earlier. For more details see the structure_tensor_analysis_v1 code. Below we will show have to combine the aggregated metrics and display them as histograms.
Another option is to save all the metrics to disk (using memory-maps as described earlier) with reasonable precision and perform the statistics directly on the full orientation volumes. This approach would be similar to what's done in the StructureTensorFiberAnalysisDemo notebook, except there data is kept in memory instead of on the disk. If you use this approach you may simply ignore the aggregated data returned by structure_tensor_analysis_v1, but it will require you to use a significant amount of disk space to store all the orientation and segmentation data. If the volumes are very big, working with them may also be slow and memory intensive, even if you use memmap to access the data.
Histograms¶
The histograms and class fractions shown here are similar to the ones in StructureTensorFiberAnalysisDemo and StructureTensorFiberAnalysisAdvancedDemo, except for one key difference. These are based on the already aggregated segmentation and orientation metrics from the blocks. This allows us to create the histograms faster and save disk space, as we don't have to save the large volumes to disk. However, there are also limitations to this approach, and if we want to save the ST results for later, saving the results per voxel (as done in StructureTensorFiberAnalysisAdvancedDemo) is more versitile than saving the block aggregated data.
The results list contains the binned metrics for each separate block, so before we can shown these, we have to combine them. First we will filter out elements with no data (i.e. blocks with no fibers) and convert the list to a NumPy array.
res = np.asarray([r[1] for r in results if r[1] is not None], dtype=np.object)
We start by reading out the number of fiber voxels for each block. Then we change NaN values for $\eta_O$ to 0 display the properly weighted mean for each fiber class.
fiber_counts = res[:, :, 0].astype(np.float)
eta_os = res[:, :, 1].astype(np.float)
eta_os[np.isnan(eta_os)] = 0
eta_o = np.average(eta_os, weights=fiber_counts, axis=0)
for i, e in enumerate(eta_o):
c = 'All' if i == 0 else class_names[i - 1]
print(f'eta_o ({c}): {round(e, 6)}')
Summing up the class counts across all blocks allows us to calculate the class fractions for the full volume
fiber_classes = np.sum(res[:, :, 0], axis=0)
fiber_fractions = fiber_classes[1:] / fiber_classes[0]
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()
Similarly we can sum of the binned $\theta$ values across blocks. We also calculate the approx. median value based on the binned data. Note that we are using the bar function rather than the hist function, as our data is already binned.
thetas = np.sum(res[:, :, 4], axis=0)
def hist(i, ax):
xs = np.arange(len(thetas[i])) / bins_per_degree
ys = thetas[i] / (np.sum(thetas[i]) / bins_per_degree)
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 = int(round(xs[np.argmax(np.cumsum(ys) / 4 > 0.5)]))
ax.axvline(xs[np.argmax(np.cumsum(ys) / 4 > 0.5)], c='b', ls='-', label=f'$Med={m}\degree$')
plt.legend()
ax.set_title(title)
ax.bar(xs, ys, align='edge', width=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()
if save_figures:
# Saving figures.
for i in range(len(class_names) + 1):
fig, ax = plt.subplots(1, 1, figsize=(5, 3))
hist(i, ax)
ax.set_title(None)
plt.savefig(os.path.join(fig_path, f'theta_{i}.pdf'), bbox_inches='tight')
plt.close()
To calculate mean and standard deviation properly for the aggregated data, we create two function. We need these to calculate the mean and standard deviation for the complete volume, based on the means and standard deviations for each block, which is what we have available in the res array.
def calc_mean(means, fiber_counts):
return np.average(means.astype(np.float), weights=fiber_counts, axis=0)
def calc_std(means, stds, fiber_counts):
mean = calc_mean(means, fiber_counts)
return np.sqrt((np.sum(fiber_counts * stds**2, axis=0) + np.sum(fiber_counts * (means - mean)**2, axis=0))
/ np.sum(fiber_counts, axis=0))
Using these function, we calculate mean and standard deviation for orientations in the XY-plane.
in_xy_angle_means = calc_mean(res[:, :, 5].astype(np.float), fiber_counts)
in_xy_angle_stds = calc_std(res[:, :, 5].astype(np.float), res[:, :, 6].astype(np.float), fiber_counts)
in_xy_angle = np.sum(res[:, :, 7], axis=0)
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 - 30, class_value + 30])
ax.axvline(class_value, c='k', ls='--', label=f'${class_value}\degree$')
xs = (np.arange(len(in_xy_angle[i])) - len(in_xy_angle[i]) // 2) / bins_per_degree
ys = in_xy_angle[i] / (np.sum(in_xy_angle[i]) / bins_per_degree)
if class_value == 90:
mask = np.abs(xs - class_value) > 90
xs[mask] = 180 + xs[mask]
ax.set_xticks(range(class_value - 30, class_value + 31, 10))
ax.set_xticklabels([60, 70, 80, '$\pm90$', -80, -70, -60])
ax.set_title(title)
ax.bar(xs, ys, align='edge', width=1 / bins_per_degree)
mean = in_xy_angle_means[i]
std = in_xy_angle_stds[i]
ax.set_xlabel('Angle (deg.)')
ax.set_ylabel('Fraction')
if i != 0 and class_value != 90:
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('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()
if save_figures:
# Saving figures.
for i in range(len(class_names) + 1):
fig, ax = plt.subplots(1, 1, figsize=(5, 3))
hist(i, ax)
ax.set_title(None)
plt.savefig(os.path.join(fig_path, f'xy-plane_{i}.pdf'), bbox_inches='tight')
plt.close()
And out-of-plane angle from the XY-plane.
out_xy_angle_means = calc_mean(res[:, :, 8].astype(np.float), fiber_counts)
out_xy_angle_stds = calc_std(res[:, :, 8].astype(np.float), res[:, :, 9].astype(np.float), fiber_counts)
out_xy_angle = np.sum(res[:, :, 10], axis=0)
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)
xs = (np.arange(len(out_xy_angle[i])) - len(out_xy_angle[i]) // 2) / bins_per_degree
ys = out_xy_angle[i] / (np.sum(out_xy_angle[i]) / bins_per_degree)
mean = out_xy_angle_means[i]
std = out_xy_angle_stds[i]
ax.set_xlim([-10, 10])
ax.bar(xs, ys, align='edge', width=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()
if save_figures:
# Saving figures.
for i in range(len(class_names) + 1):
fig, ax = plt.subplots(1, 1, figsize=(5, 3))
hist(i, ax)
ax.set_title(None)
plt.savefig(os.path.join(fig_path, f'xz-plane_{i}.pdf'), bbox_inches='tight')
plt.close()
Slices¶
As we chose to save the eigenvectors for each voxel in the data volume to the disk using the vec memory map, we can actually access all the eigenvectors along with the original data if we like. We also have orientation metrics avilable through the theta, in_xy_angle, out_xy_angle memory-maps, but they were saved with low precision intended only for visualization in other applications and we will not be using these data in this notebook.
Let's grab som data and their ST eigenvectors. The vectors were saved as float16 to save space, but before we can plot them we have to change their type. Generally, CPUs do not support float16 natively, so using this data type for anything but storing large data sets is generally not advisable.
Note: We've chosen some specific slices from our data, such as Z-index 40 and 80. If you Z-dimension size is 80 or less this will give an error. You can easily change which slices will be shown in the figures below, by changing/adding/removing the data_slices and vec_slices indices.
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[:, :, 40].copy(),
data[:, :, 80].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[:, :, :, 40].astype(np.float32),
vec[:, :, :, 80].astype(np.float32),
]
Let's can have a look at the data and the first element of the eigenvectors.
fig, axs = plt.subplots(1, len(data_slices), figsize=(20, 20))
for img, ax in zip(data_slices, axs):
ax.imshow(img, cmap='gray')
plt.show()
fig, axs = plt.subplots(1, len(vec_slices), figsize=(20, 20))
for img, ax in zip(vec_slices, axs):
ax.imshow(img[0])
plt.show()
Slices with angle information¶
Just as we did in StructureTensorFiberAnalysisDemo, we can use the calculate_angles function to calculate the segmentation and orientation metrics for the slices we just read from the volumes.
First we define some plotting functions, then we calculate vec_class, eta_os, theta, out_xy_angle and in_xy_angle for each of the slices and plot them. All the figures will be saved to the figures folder we created in the beginning. It easy to create figures for other slices simply be changing which slices are extracted from the volumes.
def add_scalebar(ax, length=100, scale=1, unit='μm', text=None):
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)
def savefig(ax, i):
title = ax.get_title()
ax.set_title(None)
plt.savefig(os.path.join(fig_path, f'{title}_{i}.svg'), bbox_inches='tight', pad_inches=0, dpi=300)
plt.savefig(os.path.join(fig_path, f'{title}_{i}.pdf'), bbox_inches='tight', pad_inches=0, dpi=300)
ax.set_title(title)
NOTE: You can easily change which figures are made by commenting out the calls to fig_with_colorbar in the show_metrics function.
def fig_with_colorbar(i, d, o, title, alpha=0.5, cmap=None, vmin=None, vmax=None):
"""Creates a figure with data, overlay and color bar."""
o = np.rot90(o, k=-1)
d = np.rot90(d, k=-1)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='3%', 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, length=2000, scale=voxel_size, text='2 mm')
if save_figures:
savefig(ax, i)
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(i, d, vec_class, 'Class', alpha=0.7)
fig_with_colorbar(i, d, eta_os, '$\eta_O$(X)', alpha=0.5, vmin=0, vmax=1)
fig_with_colorbar(i, d, theta, 'Angle from X (0-90 deg.)', alpha=0.5, vmin=0, vmax=90)
fig_with_colorbar(i, d, in_xy_angle, 'In-plane orientation (-90-90 deg.)', cmap='bwr', alpha=0.5, vmin=-90, vmax=90)
fig_with_colorbar(i, d, out_xy_angle, 'Out-of-plane orientation (-90-90 deg.)', cmap='bwr', alpha=0.5, vmin=-90, vmax=90)
fig_with_colorbar(i, d, theta, 'Angle from X (0-10 deg.)', alpha=0.5, vmin=0, vmax=10)
fig_with_colorbar(i, d, in_xy_angle, 'In-plane orientation (-10-10 deg.)', cmap='bwr', alpha=0.5, vmin=-10, vmax=10)
fig_with_colorbar(i, 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, fiber_threshold=fiber_threshold)
