This algorithm uses skewness in fourier transform to calculate direction of fibers in image. It does so by fitting a line to points above a given threshold in the fourier spectrum.
In [1]:
%pylab inline
from skimage.io import imread
from skimage.external import tifffile #only in v0.11
Turn of matplotlib interpolation, set gray colormap and figure size
In [2]:
plt.rcParams['image.interpolation'] = 'none'
plt.rcParams['image.cmap'] = 'gray'
figsize(4,4)
Open a console to do temporary work
In [186]:
%qtconsole
Create a sample image
In [3]:
size = 256
img = np.zeros((size,size), dtype=np.uint8)
t = linspace(start=0, stop=50*pi, endpoint=False, num=size)
x,y = meshgrid(t, t)
img[:,:] = 127 + 127*sin(x+0.2*y)
imshow(img);
Do fourier transform and find a suitable threshold
In [4]:
F = log(abs(fftshift(fft2(img))))
thresholded_F = np.copy(F)
threshold = F.max() * 0.25
thresholded_F[F<threshold] = 0
f, axs = subplots(ncols=2, figsize=(8,4))
axs[0].imshow(F)
axs[1].imshow(thresholded_F);
Find threshold based on number of pixels (highest 1‰ of them)
In [5]:
from skimage.exposure import cumulative_distribution
cdf = cumulative_distribution(F)
plot(cdf[1], cdf[0])
limit = np.where(cdf[0] > 0.999)[0].min()
threshold = cdf[1][limit]
plot([threshold, threshold], [0,1]) # vertical line
threshold
Out[5]:
In [6]:
thresholded_F = np.copy(F)
thresholded_F[F<threshold] = 0
f, axs = subplots(ncols=2, figsize=(8,4))
axs[0].imshow(F)
axs[1].imshow(thresholded_F);
Find a line which fits the data
In [7]:
# points
xx,yy = np.where(thresholded_F)
# y = mx + c -> matrix equation: y = Ap where A = [[x 1]]
# and p = [[m], [c]]
A = np.vstack([xx, np.ones_like(xx)]).T
p = np.linalg.lstsq(A, yy)[0]
# polynomial generator
poly = np.poly1d(p)
# plot it
line = [poly(0), poly(F.shape[1])]
plot([0,F.shape[1]], line)
xlim(0, F.shape[1])
ylim(F.shape[0], 0);
In [11]:
# assume the line is passing through center
angle = arctan(p[0]) / pi * 180
print(angle)
arctan(-0.2/1) / pi * 180
Out[11]:
==> we should have angle $\approx$ +80$^\circ$
Lets make degrees go counterclockwise instead (it's hard to "think around" convention):
$\theta = tan^{-1}(\frac{-y}{x})$
And convert to 0-180 range:
$\theta^` = \theta \% 180$
In [12]:
# turn left 90 degrees
angle = arctan(-p[0]) / pi * 180 + 90
# range 0-180
angle = angle % 180
angle
Out[12]:
Make a function of it all
In [13]:
def ft_line_fit_angle(img, threshold=0.999, debug=False):
"""
Calculate preferred orientation in image with a line fit in FT.
Parameters
----------
threshold : float
Percentage of pixels to include in FT for calculating
threshold. 0.999 * 512**2 = 262 pixels
Returns
-------
float
Angle
"""
# FT power spectrum
F = np.abs(fftshift(fft2(img)))
# do not calculate log(0)
F[F!=0], F[F==0] = log(F[F!=0]), log(F[F!=0].min())
# threshold
cdf = cumulative_distribution(F)
limit = np.where(cdf[0] > threshold)[0].min()
threshold_value = cdf[1][limit]
F = F > threshold_value
# points
y,x = np.where(F)
# special case (trivial solution)
if x.min() == x.max():
# we have a vertical line
angle = 90
b = [0, 1]
else:
# y = mx + c -> matrix equation: Ab = y where A = [[x 1]]
# and b = [[m], [c]]
A = np.vstack([x, np.ones_like(x)]).T
b = np.linalg.lstsq(A, y)[0]
# calculate angle (assume line goes through center)
angle = (90 - arctan(b[0]) / pi * 180) % 180
# show image, FT and fit
if debug:
f, ax = subplots(ncols=2, figsize=(8,4))
ax[0].imshow(img)
ax[1].imshow(F)
# add calculated line
# polynomial generator
p = np.poly1d(b)
height, width = img.shape
if angle != 90:
line = ([0, width], [p(0), p(width)])
else:
line = ([width//2, width//2], [0,height])
ax[1].plot(*line)
ax[1].set_title('computed angle: {:3.0f}'.format(angle))
ax[1].set_xlim(0,width)
ax[1].set_ylim(height,0)
return angle
Test it with another orientation of "fibers"
In [14]:
angle = 4
dx = sin(angle/180*pi)
dy = cos(angle/180*pi)
img = 127 + 127*sin(y*dy+x*dx)
ft_line_fit_angle(img, debug=True)
Out[14]:
A range of angles
In [484]:
for angle in range(0,180, 10):
dx = sin(angle/180*pi)
dy = cos(angle/180*pi)
img = 127 + 127*sin(y*dy+x*dx)
ft_line_fit_angle(img, debug=True)
Lets try with an image
In [539]:
from itertools import product
img = imread('mp.png')
bs = 100 # block size
yy,xx = img.shape
blocks_y, blocks_x = yy//block_size, xx//block_size
count = 0
for j,i in product(range(blocks_y), range(blocks_x)):
x = j*bs
y = i*bs
temp_img = img[y:y+bs, x:x+bs]
if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
continue
ot = threshold_otsu(temp_img)
if (ot < 4 # threshold above noise
or (temp_img > ot).sum() < 20): # at least 20 pixels
continue
count += 1
if count > 10:
break
#figure(figsize=(2,2))
#imshow(temp_img)
#title('{:2.2f} {:2.2f}'.format(temp_img.mean(), threshold_otsu(temp_img)))
ft_line_fit_angle(temp_img, threshold=0.9, debug=True)
Iterate over blocks of the image, and make a histogram of it all
In [14]:
from itertools import product
from skimage.filters import threshold_otsu
img = imread('mp.png')
bs = 100 # block size
blocks_y = img.shape[0]//bs
blocks_x = img.shape[1]//bs
# spread spare pixels
bsy = img.shape[0] // blocks_y
bsx = img.shape[1] // blocks_x
h = np.zeros(180) # histogram
for j,i in product(range(blocks_y), range(blocks_x)):
x = j*bsx
y = i*bsy
temp_img = img[y:y+bsy, x:x+bsx]
if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
# small image
continue
ot = threshold_otsu(temp_img)
if (ot < 4): # threshold below noise
#or (temp_img > ot).sum() < 20): # at least 20 pixels
continue
angle = ft_line_fit_angle(temp_img, threshold=0.9)
angle = round(angle)
h[angle] += 1
Plot histogram
In [15]:
figsize(8,4)
plot(h);
Basic algorithm seems ok, but we see that angles are distributed around 90$^\circ$, which seems odd. Find cause and optimize the algorithm in optimize ft line fit.ipynb.