A small example demonstrating the use of 3D structure tensor for visualizing and clustering dominant orientation.
%% Cell type:code id: tags:
``` python
%matplotlib notebook
import numpy as np
import scipy.io
import st3d
```
%% Cell type:markdown id: tags:
## Reading the data
The data is a small cube from a volumetric image og fibre composite. The cube contains five bundles in layers: UD fibre (0 deg), crossing fibre (45 deg), 90 deg fibre, -45 deg bundle, UD bundle (0 deg).
## Computing the structure tensor and the dominant orientation
Computation of structure tensor requires only two parameters: the noise scale sigma and the integration scale rho. Parameter sigma controls smothing while computing gradientes, and structures smaller than sigma will be removed by smoothing. Parameter rho gives the size over the neighborhood in which the orientation is to be analysed for every volume voxel. Larger rho will result in a smoother orientation field.
Structure tensor is a 3x3 matrix, but as it is symmetrical we only carry values of 6 elements: $s_{xx}$, $s_{yy}$, $s_{zz}$, $s_{xy}$, $s_{xz}$ and $s_{yz}$.
Eigenvalues (val) carry the information about the degree of anisotropy - this is not used or visualized here. Eigenvectors (vec) carry the orientation information, as $x$, $y$, and $z$ component of the orientation vector.
%% Cell type:code id: tags:
``` python
sigma = 0.5; # noise scale
rho = 2; # integration scale
S = st3d.structure_tensor(volume, sigma, rho)
val,vec = st3d.eig_special(S)
print(f'The volume has a shape {volume.shape}, i.e. {volume.size} voxels.')
print(f'Structure tensor information is carried in a {S.shape} array.')
print(f'Orientation information is carried in a {vec.shape} array.')
```
%% Output
The volume has a shape (150, 150, 150), i.e. 3375000 voxels.
Structure tensor information is carried in a (6, 3375000) array.
Orientation information is carried in a (3, 3375000) array.
%% Cell type:markdown id: tags:
## Visualizing the dominant orientation
Here we show only dominant orientation, ignoring shape measures.
Since slicng is in $z$ direction, the arrows show $x$ and $y$ component of orientation vectors, on top of every slice.
To interactively inspect the volume use arrow keys.
When visualizing the orientation as color, I choose fan coloring which uses hsv colors for orientatins in a certain plane, here a $xy$ plane,, and gray color for the orientations orthogonal to this plane. This coloring is convenient since all fibre bundles lay in $xy$ plane.
A tensor-vector distance can be computed for every voxel. Here I compute the distance from the z-direction. In the rotated volume, a 90 deg. bundle (the middle layer) is aligned with z-direction, and has the smallest distance.
%% Cell type:code id: tags:
``` python
u = np.array([0,0,1])
dist = st3d.tensor_vector_distance(S,u)
st3d.show_vol(dist.reshape(volume.shape))
```
%% Output
%% Cell type:markdown id: tags:
## K-means like clustering of structure tensor using tensor-vector distance
Tensor-vector distance allows clustering tensors using an approach similar to k-means. Here, every cluster is characterized by an orientation vector, and consists of the tensors which have a small distance to this orientation. The advantage of the approach is that it operates on tensors, and does not require their eigendecomposition. Only computation of cluster orientation requires eigendecomposition.
%% Cell type:code id: tags:
``` python
t = np.sqrt(1/2)
u_clusters = np.array([[1,0,0],[0,0,1],[t,0,t],[-t,0,t]]).T # initalizing orientations close to desired solution
dist = st3d.tensor_vector_distance(S,u_clusters)
assignment = np.argmin(dist,axis = 1)
S_clusters = np.zeros((6,u_clusters.shape[1]))
for r in range(10): # iterations of k-means
for i in range(u_clusters.shape[1]): # collecting ST for all voxels in a cluster
