Newer
Older
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 29 11:30:17 2019
@author: vand@dtu.dk
"""
import numpy as np
import scipy.io
import scipy.ndimage
import matplotlib.pyplot as plt
import st3d
# reading in the data
volume = scipy.io.loadmat('example_data_3D/multi_cube.mat')['vol']
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
# computing structure tensor and orientations
sigma = 0.5;
rho = 2;
S = st3d.structure_tensor(volume, sigma, rho)
val,vec = st3d.eig_special(S)
# visualization - use arrow keys to inspect the volume
st3d.show_vol_flow(volume, vec[0:2].reshape((2,)+volume.shape), s=10, double_arrow = True)
st3d.show_vol_orientation(volume, vec, coloring = st3d.fan_coloring)
#%% AND FOR A ROTATED SAMPLE...
volume = np.transpose(volume, (1, 0, 2))
S = st3d.structure_tensor(volume, sigma, rho)
val,vec = st3d.eig_special(S)
coloring_z = lambda v : st3d.fan_coloring(v[[2, 0, 1]]) # rotating coloring
st3d.show_vol_orientation(volume, vec, coloring = coloring_z)
#%% DISTANCE FROM A DIRECTION
u = np.array([0,0,1])
dist = st3d.tensor_vector_distance(S,u)
st3d.show_vol(dist.reshape(volume.shape))
# #%% FLOW MAKES SENSE ONLY WITH ROUGHLY UD FIBERS
# fxy = solve_flow_2D(S).reshape((2,)+volume.shape)
#%% K-MEANS LIKE CLUSTERING, INITILIZED BY 4 DIRECTIONS
t = np.sqrt(1/2)
u_clusters = np.array([[1,0,0],[0,0,1],[t,0,t],[-t,0,t]]).T # initalizing 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
S_clusters[:,i] = np.mean(S[:,assignment==i],axis=1)
val,vec = st3d.eig_special(S_clusters) # estimating new cluster orientation
print(f'Iter {r}: moved for {np.sum(u_clusters-vec)}')
u_clusters = vec
dist = st3d.tensor_vector_distance(S,u_clusters)
assignment = np.argmin(dist,axis = 1)
st3d.show_vol(assignment.reshape(volume.shape), cmap='jet')