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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Anders Bjorholm Dahl
abda@dtu.dk
2020
"""
import numpy as np
import skimage.feature
import cv2
# Gaussian derivatives
def getGaussDerivative(t):
'''
Computes kernels of Gaussian and its derivatives.
Parameters
----------
t : float
Vairance - t.
Returns
-------
g : numpy array
Gaussian.
dg : numpy array
First order derivative of Gaussian.
ddg : numpy array
Second order derivative of Gaussian
dddg : numpy array
Third order derivative of Gaussian.
'''
kSize = 5
s = np.sqrt(t)
x = np.arange(int(-np.ceil(s*kSize)), int(np.ceil(s*kSize))+1)
x = np.reshape(x,(-1,1))
g = np.exp(-x**2/(2*t))
g = g/np.sum(g)
dg = -x/t*g
ddg = -g/t - x/t*dg
dddg = -2*dg/t - x/t*ddg
return g, dg, ddg, dddg
# Show circles
def getCircles(coord, scale):
'''
Comptue circle coordinages
Parameters
----------
coord : numpy array
2D array of coordinates.
scale : numpy array
scale of individual blob (t).
Returns
-------
circ_x : numpy array
x coordinates of circle. Each column is one circle.
circ_y : numpy array
y coordinates of circle. Each column is one circle.
'''
theta = np.arange(0, 2*np.pi, step=np.pi/100)
theta = np.append(theta, 0)
circ = np.array((np.cos(theta),np.sin(theta)))
n = coord.shape[0]
m = circ.shape[1]
circ_y = np.sqrt(2*scale)*circ[[0],:].T*np.ones((1,n)) + np.ones((m,1))*coord[:,[0]].T
circ_x = np.sqrt(2*scale)*circ[[1],:].T*np.ones((1,n)) + np.ones((m,1))*coord[:,[1]].T
return circ_x, circ_y
# Function for detecting fibers
def detectFibers(im, diameterLimit, stepSize, tCenter, thresMagnitude):
'''
Detects fibers in images by finding maxima of Gaussian smoothed image
Parameters
----------
im : numpy array
Image.
diameterLimit : numpy array
2 x 1 vector of limits of diameters of the fibers (in pixels).
stepSize : float
step size in pixels.
tCenter : float
Scale of the Gaussian for center detection.
thresMagnitude : float
Threshold on blob magnitude.
Returns
-------
coord : TYPE
DESCRIPTION.
scale : TYPE
DESCRIPTION.
'''
radiusLimit = diameterLimit/2
radiusSteps = np.arange(radiusLimit[0], radiusLimit[1]+0.1, stepSize)
tStep = radiusSteps**2/np.sqrt(2)
r,c = im.shape
n = tStep.shape[0]
L_blob_vol = np.zeros((r,c,n))
for i in range(0,n):
g, dg, ddg, dddg = getGaussDerivative(tStep[i])
L_blob_vol[:,:,i] = tStep[i]*(cv2.filter2D(cv2.filter2D(im,-1,g),-1,ddg.T) +
cv2.filter2D(cv2.filter2D(im,-1,ddg),-1,g.T))
# Detect fibre centers
g, dg, ddg, dddg = getGaussDerivative(tCenter)
Lg = cv2.filter2D(cv2.filter2D(im, -1, g), -1, g.T)
coord = skimage.feature.peak_local_max(Lg, threshold_abs = thresMagnitude)
# Find coordinates and size (scale) of fibres
magnitudeIm = np.min(L_blob_vol, axis = 2)
scaleIm = np.argmin(L_blob_vol, axis = 2)
scales = scaleIm[coord[:,0], coord[:,1]]
magnitudes = -magnitudeIm[coord[:,0], coord[:,1]]
idx = np.where(magnitudes > thresMagnitude)
coord = coord[idx[0],:]
scale = tStep[scales[idx[0]]]
return coord, scale