In [109]:
%pylab inline
import skimage.io as io
plt.rcParams['image.cmap'] = 'gray_r'
plt.rcParams['image.interpolation'] = 'none'
In [13]:
img = io.imread('1_66500.png')
img.shape
Out[13]:
In [17]:
files1 = !ls 1_*.png
files2 = !ls 2_*.png
len(files1), len(files2)
Out[17]:
In [18]:
for n, filename in enumerate(files1):
if n%2 == 0:
fig, axs = subplots(nrows=1, ncols=2)
img = io.imread(filename);
axs[n%2].imshow(img[:,:,0])
axs[n%2].set_title('class 1 - ' + filename)
In [19]:
for n, filename in enumerate(files2):
if n%2 == 0:
fig, axs = subplots(nrows=1, ncols=2)
img = io.imread(filename);
axs[n%2].imshow(img[:,:,0])
axs[n%2].set_title('class 2 - ' + filename)
Wrong channel for 2_88* ?
In [20]:
img = io.imread('2_88100.png')
img.shape
Out[20]:
In [22]:
for i in range(3):
figure()
imshow(img[:,:,i])
In [29]:
%qtconsole
In [25]:
!mkdir original
!mv *.png original
In [41]:
files1 = !ls original/1_*png
files2 = !ls original/2_6*png
files1.extend(files2)
for f in files1:
img = io.imread(f)[:,:,0]
io.imsave(f.replace('original/', ''), img);
In [42]:
files2 = !ls original/2_8*png
for f in files2:
img = io.imread(f)[:,:,1]
io.imsave(f.replace('original/', ''), img);
In [4]:
from skimage.exposure import cumulative_distribution
from scipy.stats import linregress
from skimage import morphology
def angle_ft_line_fit(img, debug=False):
"""
Calculate preferred orientation in image with a line fit
in FT.
Parameters
----------
img : ndarray
Image to do FT and find line fit in FT.
debug : bool
If function should show a image for debugging.
Returns
-------
(angle, r) : tuple of floats
Angle of line fit and r value of line fit.
"""
# 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())
# only include top value pixels in FT
number_of_line_pixels = img.shape[0]//3 # 33% of height
threshold = 1 - number_of_line_pixels / img.size
# use cdf to find maximum pixels
cdf = cumulative_distribution(F)
limit = np.where(cdf[0] > threshold)[0].min()
threshold_value = cdf[1][limit]
F = F > threshold_value #binary image
# create convex hull, to boost correlation of lines vs blobs
F = morphology.convex_hull_image(F)
# points
y,x = np.where(F)
# if dy > dx -> swap x,y in line fit for accuracy
dx = abs(x.max()-x.min())
dy = abs(y.max()-y.min())
# solve x=my+c by least-squares regression
if dx < dy:
m,c,r,pv,err = linregress(y,x)
b = (1/m, -c/m)
# calculate angle - assume line goes through center
angle = (90 - arctan(b[0]) / pi * 180) % 180
# solve y=mx+c by least-squares regression
else:
line_span = dx
m,c,r,pv,err = linregress(x,y)
b = (m,c)
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 != 0:
line = ([0, width], [p(0), p(width)])
else:
line = ([width//2, width//2], [0,height])
ax[1].plot(*line)
ax[1].set_title('ang: {:3.0f} r:{:0.2} err:{:0.2}'
.format(angle,r,err))
ax[1].set_xlim(0,width)
ax[1].set_ylim(height,0)
return (angle, r)
In [68]:
from itertools import product
from skimage.filter import threshold_otsu
img = io.imread('1_66500.png')
iy,ix = img.shape
by,bx = (5,5) # number of blocks in x,y direction
bs = img.shape[0] // by
bsy, bsx = iy//by, ix//bx # block size
count = 0
f, axs = subplots(nrows=3, ncols=4, figsize=(10,8))
for j,i in product(range(by), range(bx)):
x,y = j*bs, i*bs
temp_img = img[y:y+bs, x:x+bs]
mean = temp_img.mean()
if mean <= 5:
continue
ot = threshold_otsu(temp_img)
if ot < 1:
continue
if count >= 12:
break
ax = axs[count//4, count%4]
ax.imshow(temp_img)
ax.set_title('m: {:1.1f}, t:{:1.1f}'.format(mean, ot))
count += 1
if count == 6: # pick out image for manual debug
ii = np.copy(temp_img)
angle_ft_line_fit(temp_img, debug=True)
In [70]:
img = io.imread('2_66500.png')
iy,ix = img.shape
by,bx = (5,5) # number of blocks in x,y direction
bs = img.shape[0] // by
bsy, bsx = iy//by, ix//bx # block size
count = 0
f, axs = subplots(nrows=3, ncols=4, figsize=(10,8))
for j,i in product(range(by), range(bx)):
x,y = j*bs, i*bs
temp_img = img[y:y+bs, x:x+bs]
if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
continue
mean = temp_img.mean()
if mean <= 5:
continue
ot = threshold_otsu(temp_img)
if ot < 1:
continue
if count >= 12:
break
ax = axs[count//4, count%4]
ax.imshow(temp_img)
ax.set_title('m: {:1.1f}, t:{:1.1f}'.format(mean, ot))
count += 1
if count == 6: # pick out image for manual debug
ii = np.copy(temp_img)
angle_ft_line_fit(temp_img, debug=True)
In [26]:
from skimage.exposure import cumulative_distribution
from scipy.stats import linregress
from skimage import morphology
%pylab inline
def angle_ft_line_fit(img, debug=False):
"""
Calculate preferred orientation in image with a line fit
in FT.
Parameters
----------
img : ndarray
Image to do FT and find line fit in FT.
debug : bool
If function should show a image for debugging.
Returns
-------
(angle, r) : tuple of floats
Angle of line fit and r value of line fit.
"""
# FT power spectrum
F = np.abs(fftshift(fft2(img)))
# do not calculate log(0)
F = np.log(F.clip(1))
# only include top value pixels in FT
number_of_line_pixels = img.shape[0]//3 # 33% of height
threshold = 1 - number_of_line_pixels / img.size
# use cdf to find maximum pixels
cdf = cumulative_distribution(F)
limit = np.where(cdf[0] > threshold)[0].min()
threshold_value = cdf[1][limit]
F = F > threshold_value #binary image
# create convex hull, to boost correlation of lines vs blobs
F = morphology.convex_hull_image(F)
# points
x,y = np.where(F)
# if dy > dx -> swap x,y in line fit for accuracy
dx = abs(x.max()-x.min())
dy = abs(y.max()-y.min())
# solve x=my+c by least-squares regression
if dx < dy:
m,c,r,pv,err = linregress(y,x)
b = (1/m, -c/m)
# calculate angle - assume line goes through center
angle = (90 - arctan(b[0]) / pi * 180) % 180
# solve y=mx+c by least-squares regression
else:
m,c,r,pv,err = linregress(x,y)
b = (m,c)
angle = (90 - arctan(b[0]) / pi * 180) % 180
return (angle, r)
def feature1(img_filename, r_above=0.8):
"""
Returns number of subimages with r value above 0.8 for
angle_ft_line_fit()
"""
from skimage import io
from itertools import product
img = io.imread(img_filename);
iy,ix = img.shape
by,bx = (5,5) # number of blocks in x,y direction
bsy,bsx = iy//by, ix//bx # block size
n = 0
for i,j in product(range(by), range(bx)):
y,x = i*bsy, j*bsx
temp_img = img[y:y+bsy, x:x+bsx]
# discard empty subimages
mean = temp_img.mean()
if mean <= 5:
continue
a, r = angle_ft_line_fit(temp_img)
if r > r_above:
n+=1
return n
In [27]:
%time feature1('1_66500.png'), feature1('1_88100.png'), feature1('2_66500.png')
Out[27]:
In [140]:
files = !ls [0-2]*png
y = []
for f in files:
y.append(feature1(f))
plot(y)
Out[140]:
Might also work for whole image
In [30]:
from skimage import io
img = io.imread('fibers/1_66500.png')
angle_ft_line_fit(img, debug=True)
img = io.imread('fibers/1_88100.png')
angle_ft_line_fit(img, debug=True)
img = io.imread('fibers/2_66500.png')
angle_ft_line_fit(img, debug=True)
Out[30]:
In [119]:
%pylab inline
from skimage import io, filters
In [31]:
%qtconsole
In [114]:
def ft_radial_1d_signature(img):
"""
Compute radial signature of FT where values are above a threshold.
"""
# Fourier spectrum
F = np.abs(fftshift(fft2(img)))
F = np.log(F.clip(1))
# threshold
F[F<F.max()*0.65] = 0
# positions
x,y = np.where(F!=0)
# center of mass
# should be center, but its easy to calculate
cm = y.mean(), x.mean()
# make (0,0) to cm
y = y - cm[0]
x = x - cm[1]
angles = np.zeros(180)
n = 0
for point in zip(y,x):
# only upper half, as its symetric
if point[0] < 0:
continue
# special case in range 0-180
if point[0] == 0 and point[1] < 0:
continue
# center of mass
if point == (0,0):
continue
angle = floor(arctan2(*point) / pi * 180)
#angles[angle] += F[point[0]-cm[0], point[1]-cm[1]]
angles[angle] += 1
return angles, F
In [123]:
def _f(filename):
img = io.imread(filename)
h, F = ft_radial_1d_signature(img)
figure(figsize=(12,4))
subplot(131)
imshow(img)
subplot(132)
imshow(F)
subplot(133)
plot(range(180), h)
_f('1_66500.png')
_f('2_66500.png')
_f('1_88100.png')
_f('2_88100.png')
In [ ]:
png
In [ ]:
from scipy.optimize import curve_fit
In [147]:
def gaussian(x, m, sigma, C):
return C * e**(-0.5*(x-m)**2/sigma**2)
x = np.linspace(0,10)
y = gaussian(x, 3, 1, 10)
plot(x,y);
Out[147]:
In [170]:
h, F = ft_radial_1d_signature(img)
fit = curve_fit(gaussian, range(180), h)[0]
print(fit)
plot(h)
plot(gaussian(range(180), *fit))
Out[170]:
In [125]:
from skimage import io, filters
from scipy.optimize import curve_fit
def feature2(filename):
"""
Compute angle distribution of thresholded ft.
Return sigma of best gaussian fit to angle distribution.
"""
def ft_radial_1d_signature(filename):
"""
Compute radial signature of FT where values are above a threshold.
"""
img = io.imread(filename)
# Fourier spectrum
F = np.abs(fftshift(fft2(img)))
F = np.log(F.clip(1))
# threshold
F[F<F.max()*0.65] = 0
# positions
x,y = np.where(F!=0)
# center of mass
# should be center, but its easy to calculate
cm = y.mean(), x.mean()
# make (0,0) to cm
y = y - cm[0]
x = x - cm[1]
angles = np.zeros(180)
n = 0
for point in zip(y,x):
# only upper half, as its symetric
if point[0] < 0:
continue
# special case in range 0-180
if point[0] == 0 and point[1] < 0:
continue
# center of mass
if point == (0,0):
continue
angle = floor(arctan2(*point) / pi * 180)
#angles[angle] += F[point[0]-cm[0], point[1]-cm[1]]
angles[angle] += 1
return angles
def gaussian(x, m, sigma, C):
return C * e**(-0.5*(x-m)**2/sigma**2)
# get angle histogram from ft
h = ft_radial_1d_signature(filename)
# compute best gaussian fit
try:
m, sigma, C = curve_fit(gaussian, range(180), h)[0]
except:
# no fit
return 180
if sigma > 180:
sigma = 180
return sigma
In [126]:
%time feature2('1_66500.png'), feature2('1_88100.png'), feature2('2_66500.png')
Out[126]:
In [127]:
files = !ls [0-2]*.png
y = []
for f in files:
y.append(feature2(f))
plot(y)
Out[127]:
In [187]:
from skimage import filters
In [251]:
img = io.imread('1_88100.png')
imshow(img);
Separate weak and strong signal
In [197]:
threshold = 10
timg = np.copy(img)
timg[img<threshold] = 0
timg[img>=threshold] = 1
imshow(timg);
Only keep large features
In [198]:
from skimage import morphology
In [249]:
for i in range(3,15):
oimg = morphology.binary_opening(timg, morphology.disk(i))
figure()
imshow(oimg)
colorbar();
In [253]:
for i in range(3,15):
oimg = morphology.opening(img, morphology.disk(i))
figure()
imshow(oimg)
colorbar();
In [256]:
img = io.imread('1_66500.png')
for i in range(3,15):
oimg = morphology.binary_opening(img>10, morphology.disk(i))
figure()
imshow(oimg)
colorbar();
In [257]:
img = io.imread('2_66500.png')
for i in range(3,15):
oimg = morphology.binary_opening(img>10, morphology.disk(i))
figure()
imshow(oimg)
colorbar();
In [38]:
from skimage import io, morphology, measure, filters
def feature3(filename):
"""
Find number of connected areas after an opening of a disc with
radius=3.
"""
img = io.imread(filename)
opening = morphology.binary_opening(img>10, morphology.disk(3))
labeled = measure.label(opening)
#figure()
#imshow(labeled)
#title(f)
#colorbar()
return labeled.max()
In [289]:
files = !ls *png
y = []
for f in files:
n = feature3(f)
y.append(n)
plot(y)
Out[289]:
In [39]:
from skimage import io, morphology, measure, filters
def feature4(filename, r=10):
"""
Find number of connected areas after an dilation of a disc with
radius=10.
"""
img = io.imread(filename)
opening = morphology.binary_dilation(img>10, morphology.disk(r))
labeled = measure.label(opening)
#figure()
#imshow(labeled)
#title(f)
#colorbar()
return labeled.max()
In [292]:
files = !ls *png
y = []
for f in files:
n = feature4(f)
y.append(n)
plot(y)
Out[292]:
In [142]:
from sklearn import datasets
from random import shuffle
In [143]:
type(datasets.load_iris())
Out[143]:
In [144]:
files = !ls [0-2]*png
# randomize
shuffle(files)
n = len(files)
feat = 4
data = np.zeros((n,feat))
target = np.zeros(n)
fibers = datasets.base.Bunch(data=data, target=target)
In [145]:
for n,f in enumerate(files):
f1 = feature1(f)
f2 = feature2(f)
f3 = feature3(f)
f4 = feature4(f)
fibers.data[n] = [f1,f2,f3,f4]
fibers.target[n] = f[0:1]
In [146]:
fibers.data[0], fibers.target[0]
Out[146]:
In [147]:
# store filename and image too
fibers.files = files
In [148]:
i = 10
fibers.files[i], fibers.data[i], fibers.target[i]
Out[148]:
save data and classifier to pickle
In [149]:
from sklearn.externals import joblib
joblib.dump(fibers, 'fibers.pkl')
Out[149]:
In [150]:
from sklearn import tree
from sklearn.externals import joblib
#from sklearn.utils import shuffle
#x,y,z = shuffle(x,y,z)
In [151]:
fibers = joblib.load('fibers.pkl')
len(fibers.files)*0.8
Out[151]:
In [152]:
clf = tree.DecisionTreeClassifier(criterion='entropy', max_depth=4)
In [153]:
clf.fit(fibers.data[0:60], fibers.target[0:60])
Out[153]:
In [154]:
joblib.dump(clf, 'clf.pkl')
Out[154]:
In [155]:
from sklearn import tree
from sklearn.externals.six import StringIO
from sklearn.externals import joblib
import pydot
In [156]:
clf = joblib.load('clf.pkl')
In [157]:
dot_data = StringIO()
tree.export_graphviz(clf, out_file=dot_data)
graph = pydot.graph_from_dot_data(dot_data.getvalue())
graph.write_png('decision_tree.png')
Out[157]:
In [158]:
from IPython.display import Image
Image('decision_tree.png')
Out[158]:
In [159]:
wrong = 0
right = 0
for i in range(60, len(fibers.files)):
predicted = clf.predict(fibers.data[i])
real = fibers.target[i]
filename = fibers.files[i]
if predicted != real:
print('WRONG! File: {}, Predicted: {}'.format(filename, predicted))
wrong += 1
else:
right += 1
print('Right: {}, Wrong: {}'.format(right, wrong))
In [162]:
wrong = []
right = 0
for i in range(len(fibers.files)):
predicted = clf.predict(fibers.data[i])
real = fibers.target[i]
filename = fibers.files[i]
if predicted != real:
print('WRONG! File: {}, Predicted: {}'.format(filename, predicted))
wrong.append(filename)
else:
right += 1
print('Right: {}, Wrong: {}'.format(right, len(wrong)))
In [163]:
w = wrong[0]
print(feature1(w), feature2(w), feature3(w), feature4(w))
Image(wrong[0])