A small example demonstrating the use of 3D structure tensor for visualizing and clustering dominant orientation.
A small example demonstrating the use of 3D structure tensor for visualizing and clustering dominant orientation.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%matplotlib notebook
%matplotlib notebook
import numpy as np
import numpy as np
import scipy.io
import scipy.io
import st3d
import st3d
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Reading the data
## 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).
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
## 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.
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}$.
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.
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:
%% Cell type:code id: tags:
``` python
``` python
sigma = 0.5; # noise scale
sigma = 0.5; # noise scale
rho = 2; # integration scale
rho = 2; # integration scale
S = st3d.structure_tensor(volume, sigma, rho)
S = st3d.structure_tensor(volume, sigma, rho)
val,vec = st3d.eig_special(S)
val,vec = st3d.eig_special(S)
print(f'The volume has a shape {volume.shape}, i.e. {volume.size} voxels.')
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'Structure tensor information is carried in a {S.shape} array.')
print(f'Orientation information is carried in a {vec.shape} array.')
print(f'Orientation information is carried in a {vec.shape} array.')
```
```
%% Output
%% Output
The volume has a shape (150, 150, 150), i.e. 3375000 voxels.
The volume has a shape (150, 150, 150), i.e. 3375000 voxels.
Structure tensor information is carried in a (6, 3375000) array.
Structure tensor information is carried in a (6, 3375000) array.
Orientation information is carried in a (3, 3375000) array.
Orientation information is carried in a (3, 3375000) array.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Visualizing the dominant orientation
## Visualizing the dominant orientation
Here we show only dominant orientation, ignoring shape measures.
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.
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.
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.
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.
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:
%% Cell type:code id: tags:
``` python
``` python
u = np.array([0,0,1])
u = np.array([0,0,1])
dist = st3d.tensor_vector_distance(S,u)
dist = st3d.tensor_vector_distance(S,u)
st3d.show_vol(dist.reshape(volume.shape))
st3d.show_vol(dist.reshape(volume.shape))
```
```
%% Output
%% Output
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## K-means like clustering of structure tensor using tensor-vector distance
## 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.
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:
%% Cell type:code id: tags:
``` python
``` python
t = np.sqrt(1/2)
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
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)
dist = st3d.tensor_vector_distance(S,u_clusters)
assignment = np.argmin(dist,axis = 1)
assignment = np.argmin(dist,axis = 1)
S_clusters = np.zeros((6,u_clusters.shape[1]))
S_clusters = np.zeros((6,u_clusters.shape[1]))
for r in range(10): # iterations of k-means
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
for i in range(u_clusters.shape[1]): # collecting ST for all voxels in a cluster
