Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#!/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.ndimage
# STRUCTURE TENSOR 2D
def structure_tensor(image, sigma, rho):
""" Structure tensor for 2D image data
Arguments:
image: a 2D array of size N = rows*columns
sigma: a noise scale, structures smaller than sigma will be
removed by smoothing
rho: an integration scale giving the size over the neighborhood in
which the orientation is to be analysed
Returns:
an array with shape (3,N) containing elements of structure tensor
s_xx, s_yy, s_xy ordered acording to image.ravel()
Author: vand@dtu.dk, 2019
"""
# computing derivatives (scipy implementation truncates filter at 4 sigma)
image = image.astype(np.float);
Ix = scipy.ndimage.gaussian_filter(image, sigma, order=[1,0], mode='nearest')
Iy = scipy.ndimage.gaussian_filter(image, sigma, order=[0,1], mode='nearest')
# integrating elements of structure tensor (scipy uses sequence of 1D)
Jxx = scipy.ndimage.gaussian_filter(Ix**2, rho, mode='nearest')
Jyy = scipy.ndimage.gaussian_filter(Iy**2, rho, mode='nearest')
Jxy = scipy.ndimage.gaussian_filter(Ix*Iy, rho, mode='nearest')
S = np.vstack((Jxx.ravel(), Jyy.ravel(), Jxy.ravel()));
return S
def eig_special(S):
""" Eigensolution for symmetric real 2-by-2 matrices
Arguments:
S: an array with shape (3,N) containing structure tensor
Returns:
val: an array with shape (2,N) containing sorted eigenvalues
vec: an array with shape (2,N) containing eigenvector corresponding
to the smallest eigenvalue (the other is orthogonal to the first)
Author: vand@dtu.dk, 2019
"""
val = 0.5*(S[0]+S[1]+np.outer(np.array([-1,1]), np.sqrt((S[0]-S[1])**2+4*S[2]**2)))
vec = np.vstack((-S[2],S[0]-val[0])) # y will be positive
aligned = S[2]==0 # dealing with diagonal matrices
vec[:,aligned] = 1-np.argsort(S[0:2,aligned], axis=0)
vec = vec/np.sqrt(vec[0]**2+vec[1]**2) # normalizing
return val, vec
def solve_flow(S):
""" Solving 1D optic flow, returns LLS optimal x for flow along y axis
( x is a solution to S[0]*x=S[2] )
Arguments:
S: an array with shape (3,N) containing 2D structure tensor
Returns:
x: an array with shape (1,N) containing x components of the flow
Author: vand@dtu.dk, 2019
"""
aligned = S[0]==0 # 0 or inf solutions
x = np.zeros((1,S.shape[1]))
x[0,~aligned] = - S[2,~aligned]/S[0,~aligned]
return x # returning shape (1,N) array for consistancy with 3D case
#% HELPING FUNCTIONS FOR VISUALIZATION
def plot_orientations(ax, dim, vec, s = 5):
""" Helping function for adding orientation-quiver to the plot.
Arguments: plot axes, image shape, orientation, arrow spacing.
"""
vx = vec[0].reshape(dim)
vy = vec[1].reshape(dim)
xmesh, ymesh = np.meshgrid(np.arange(dim[0]), np.arange(dim[1]), indexing='ij')
ax.quiver(ymesh[s//2::s,s//2::s],xmesh[s//2::s,s//2::s],vy[s//2::s,s//2::s],vx[s//2::s,s//2::s],color='r',angles='xy')
ax.quiver(ymesh[s//2::s,s//2::s],xmesh[s//2::s,s//2::s],-vy[s//2::s,s//2::s],-vx[s//2::s,s//2::s],color='r',angles='xy')
def polar_histogram(ax, distribution, cmap = 'hsv'):
""" Helping function for producing polar histogram.
Arguments: plot axes, oriantation distribution, colormap.
"""
N = distribution.size
bin_centers_full = (np.arange(2*N)+0.5)*np.pi/N # full circle (360 deg)
distribution_full = np.r_[distribution,distribution]/max(distribution) # full normalized distribution
x = np.r_[distribution_full*np.cos(bin_centers_full),0]
y = np.r_[distribution_full*np.sin(bin_centers_full),0]
triangles = np.array([(i, np.mod(i-1,2*N), 2*N) for i in range(2*N)]) # triangles[0] is symmetric over 0 degree
triangle_centers_full = (np.arange(2*N))*np.pi/N # a triangle covers area BETWEEN two bin_centers
triangle_colors = np.mod(triangle_centers_full, np.pi)/np.pi # from 0 to 1-(1/2N)
ax.tripcolor(y, x, triangles, facecolors=triangle_colors, cmap=cmap, vmin = 0.0, vmax = 1.0)
ax.set_aspect('equal')
ax.set_xlim([-1,1])
ax.set_ylim([1,-1])