%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
from scipy.misc import lena
lenabw = lena()/float(lena().max())
imshow(lenabw, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f522663dd50>
lena().max(), lena().dtype
(245, dtype('int64'))
from scipy.ndimage.filters import sobel
sx = sobel(lenabw, axis=1)
sy = sobel(lenabw, axis=0)
sxy = hypot(sx, sy)
imshow(sxy, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f52247ff290>
imshow(sxy*100.0, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f52245b8128>
from scipy.ndimage import measurements
h = measurements.histogram(lenabw, 0, 1, 256)
plot(h)
grid()
h = measurements.histogram(sxy, 0, sxy.max(), 256)
plot(h)
[<matplotlib.lines.Line2D at 0x7f5224433c10>]
plot(linspace(0, sxy.max(), len(h)), h)
grid()
subplot(131)
imshow(where(sxy > 0.2, 1, 0), cmap=cm.gray)
subplot(132)
imshow(where(sxy > 0.4, 1, 0), cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f52243662d0>
plot(linspace(0, sxy.max(), len(h)), h)
grid()
vlines((0.1, 0.2, 0.4), 0.001, h.max())
<matplotlib.collections.LineCollection at 0x7f52244f0d50>
imshow(where(sxy > 0.1, 1, 0), cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f5224234f50>
from scipy.integrate import cumtrapz
plot(cumtrapz(h), "o-")
grid()
argwhere(cumtrapz(h) > cumtrapz(h).max()/2.0)[0]
array([7])
7/256.0
0.02734375
argwhere(cumtrapz(h) > cumtrapz(h).max()*0.75)[0]/256.0
array([ 0.0703125])
To automate threshold setting, maybe you want to find the point at which the derivative of the cumulative integral is a certain value (i.e. the slope of the cumtrapz output)
Edge detection is easy by looking at the derivative of a signal:
edge = r_[zeros((200)), linspace(0,1, 100), ones((100))]
subplot(131); imshow((edge,), aspect='auto', cmap=cm.gray)
subplot(132); plot(edge)
subplot(133); plot(diff(edge))
[<matplotlib.lines.Line2D at 0x7f522407c510>]
But, what if there is noise?
noisy_edge = edge + random.random(edge.shape)*0.05
subplot(131); imshow((noisy_edge,), aspect='auto', cmap=cm.gray)
subplot(132); plot(noisy_edge)
subplot(133); plot(diff(noisy_edge))
[<matplotlib.lines.Line2D at 0x7f521f360910>]
subplot(121); plot(noisy_edge)
subplot(122); plot(convolve(noisy_edge, ones(50)/50.0)[:400])
[<matplotlib.lines.Line2D at 0x7f521f186b10>]
edge = r_[zeros((200)), linspace(0,1, 100), ones((100))]
subplot(141); imshow((edge,), aspect='auto', cmap=cm.gray)
subplot(142); plot(edge)
subplot(143); plot(diff(edge))
subplot(144); plot(diff(diff(edge)))
[<matplotlib.lines.Line2D at 0x7f521ef64c50>]
smoothed = convolve(noisy_edge, ones(100)/100.0, mode='valid')
subplot(131)
plot(smoothed)
subplot(132)
plot(diff(smoothed))
subplot(133)
plot(diff(diff(smoothed)))
ylim((-0.003, 0.003))
(-0.003, 0.003)
from scipy.ndimage import gaussian_laplace
LoG = gaussian_laplace(lenabw, 4)
imshow(LoG, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f521ebb89e0>
th = LoG > 0
subplot(132)
imshow(th, cmap=cm.gray);
from scipy.ndimage.filters import convolve as convolve2
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges = convolve2(th, kernel)
subplot(121)
imshow(MHedges, cmap=cm.gray)
colorbar()
subplot(122)
imshow(lenabw, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f521ea51200>
LoG = gaussian_laplace(lenabw, 8)
subplot(131); imshow(LoG, cmap=cm.gray)
th = LoG > 0
subplot(132); imshow(th, cmap=cm.gray)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges_lena = convolve2(th, kernel)
subplot(133); imshow(MHedges_lena, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e998850>
LoG = gaussian_laplace(lenabw, 1)
subplot(131); imshow(LoG, cmap=cm.gray)
th = LoG > 0
subplot(132); imshow(th, cmap=cm.gray)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges_lena2 = convolve2(th, kernel)
subplot(133); imshow(MHedges_lena2, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e79d250>
fruits = sum(imread('vitamy fruits.png'), axis=2)
imshow(fruits, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e616190>
LoG = gaussian_laplace(fruits, 4)
subplot(131)
imshow(LoG, cmap=cm.gray)
subplot(132)
th = LoG > 0
imshow(th, cmap=cm.gray)
subplot(133)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges = convolve2(th, kernel)
imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e4c4f90>
imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e33e290>
LoG = gaussian_laplace(fruits, 8)
subplot(131); imshow(LoG, cmap=cm.gray)
th = LoG > 0
subplot(132); imshow(th, cmap=cm.gray)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges = convolve2(th, kernel)
subplot(133); imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e274d10>
LoG = gaussian_laplace(fruits, 8)
subplot(131); imshow(LoG, cmap=cm.gray)
th = LoG > 0
subplot(132); imshow(th, cmap=cm.gray)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges = convolve2(th, kernel)
subplot(133)
imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521e078710>
imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521dee69d0>
coins = imread('coins.png')
imshow(coins, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521de20810>
LoG = gaussian_laplace(coins, 8)
subplot(131); imshow(LoG, cmap=cm.gray)
th = LoG > 0
subplot(132); imshow(th, cmap=cm.gray)
kernel = [[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]]
MHedges = convolve2(th, kernel)
subplot(133)
imshow(MHedges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521dcdef90>
from scipy.ndimage.filters import gaussian_filter
sm = gaussian_filter(fruits, 4)
sx = sobel(sm, axis=1)
sy = sobel(sm, axis=0)
sxy = hypot(sx, sy)
imshow(sxy, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521db088d0>
sx = sobel(sm, axis=1)
sy = sobel(sm, axis=0)
a = sy/sx
a = where(isnan(a), 0, a)
fruit_angle = arctan2(sy, sx)
imshow(fruit_angle)
colorbar()
fruit_angle.max(), fruit_angle.min(), fruit_angle[500, 100]
-c:3: RuntimeWarning: divide by zero encountered in true_divide -c:3: RuntimeWarning: invalid value encountered in true_divide
(3.1415927, -3.1415491, 1.1172101)
Provides extended image processing capabilities
http://scikit-image.org/docs/dev/api/skimage.html
like:
'''canny.py - Canny Edge detector
Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
Pattern Analysis and Machine Intelligence, 8:679-714, 1986
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
Pattern Analysis and Machine Intelligence, 8:679-714, 1986
William Green' Canny tutorial
http://dasl.mem.drexel.edu/alumni/bGreen/www.pages.drexel.edu/_weg22/can_tut.html
'''
import skimage.filter as filters
fruits_canny = filters.canny(fruits, 4)
imshow(fruits_canny, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521da599d0>
fruits_canny = filters.canny(fruits, 2)
imshow(fruits_canny, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521da21250>
fruits_canny = filters.canny(fruits, 2, 0.05, 0.4)
imshow(fruits_canny, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f521341c710>
fruits_canny = filters.canny(fruits, 2, 0.14, 0.18)
imshow(fruits_canny, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f5213349a90>
lena_canny = filters.canny(lenabw, 2)
subplot(131); imshow(lena_canny, cmap= cm.gray)
subplot(132); imshow(MHedges_lena, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f52132845d0>
subplot(131); imshow(coins, cmap=cm.gray)
coins_canny = filters.canny(coins, 2)
subplot(132); imshow(coins_canny, cmap= cm.gray)
coins_canny = filters.canny(coins, 4)
subplot(133); imshow(coins_canny, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f52130de750>
aerial = imread('aerial.png')
aerialbw = sum(aerial, axis=2)/3.0
imshow(aerialbw, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f520ff7fa28>
LoG = gaussian_laplace(aerialbw, 24)
subplot(131); imshow(LoG, cmap=cm.gray)
th = where(LoG > 0, 1.0, 0.0)
subplot(132); imshow(th, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520fed56d0>
sm = gaussian_filter(th , 8)
imshow(sm, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520fdb5f90>
sm = gaussian_filter(th , 4)
imshow(1 - sm, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520fcf8450>
smoothed = 1 - sm
imshow(diff(smoothed) == 0, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520fc2f790>
Then morphology filters and skeletonize.
sm = where(sm > 0.5, 0, 1)
imshow(sm, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520fb62a50>
import skimage.morphology as morphology
eroded = morphology.binary_erosion(sm, ones((20,20)))
imshow(eroded, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520e50c310>
subplot(121)
imshow(eroded[400:700, 400:700], cmap=cm.gray, interpolation='nearest')
subplot(122)
imshow(morphology.medial_axis(eroded[400:800, 400:800]), cmap=cm.gray, interpolation='nearest')
<matplotlib.image.AxesImage at 0x7f520e44d150>
#from scikits image
def hough(img, theta=None):
if img.ndim != 2:
raise ValueError('The input image must be 2-D')
if theta is None:
theta = np.linspace(-np.pi / 2, np.pi / 2, 179, endpoint=False)
# compute the vertical bins (the distances)
d = np.ceil(np.hypot(*img.shape))
nr_bins = 2 * d
bins = np.linspace(-d, d, nr_bins)
# allocate the output image
out = np.zeros((nr_bins, len(theta)), dtype=np.uint64)
# precompute the sin and cos of the angles
cos_theta = np.cos(theta)
sin_theta = np.sin(theta)
# find the indices of the non-zero values in
# the input image
y, x = np.nonzero(img)
# x and y can be large, so we can't just broadcast to 2D
# arrays as we may run out of memory. Instead we process
# one vertical slice at a time.
for i, (cT, sT) in enumerate(zip(cos_theta, sin_theta)):
# compute the base distances
distances = x * cT + y * sT
# round the distances to the nearest integer
# and shift them to a nonzero bin
shifted = np.round(distances) - bins[0]
# cast the shifted values to ints to use as indices
indices = shifted.astype(np.int)
# use bin count to accumulate the coefficients
bincount = np.bincount(indices)
# finally assign the proper values to the out array
out[:len(bincount), i] = bincount
return out, theta, bins
test_image = zeros((50, 50))
test_image[range(20),range(20)] = 1
imshow(test_image)
<matplotlib.image.AxesImage at 0x7f520d85bf10>
h, theta, d = hough(test_image)
imshow(h, cmap= cm.binary, extent=[rad2deg(theta[0]), rad2deg(theta[-1]), d[-1], d[0]])
<matplotlib.image.AxesImage at 0x7f520d787cd0>
h.shape
(142, 179)
hypot(50,50) # maximum distance measured by Hough transform
70.710678118654755
argmax(h)
12753
unravel_index(argmax(h),h.shape)
(71, 44)
print(rad2deg(theta[44]), d[71])
-45.7541899441 0.503546099291
test_image = zeros((50, 50))
test_image[range(20),range(20)] = 1
test_image[arange(24) + 15, (arange(24)*2)] = 1
imshow(test_image, cmap=cm.gray, interpolation='nearest')
ylim((0, 50))
(0, 50)
h, theta, d = hough(test_image)
imshow(h, cmap= cm.binary , extent=[rad2deg(theta[-1]), rad2deg(theta[0]), d[-1], d[0]])
unravel_index(argmax(h),h.shape)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f520d5df8c0>
h.max()
24
i_r, i_t = unravel_index(argmax(h),h.shape)
print(rad2deg(theta[i_t]), d[i_r])
-62.8491620112 -12.5886524823
img = np.zeros((100, 150), dtype=bool)
img[30, :] = 1
img[:, 65] = 1
img[35:45, 35:50] = 1
for i in range(90):
img[i, i] = 1
img += np.random.random(img.shape) > 0.95
imshow(img, cmap=cm.gray)
ylim((0,100))
(0, 100)
h, theta, d = hough(img)
imshow(h, cmap= cm.binary, extent=[rad2deg(theta[0]), rad2deg(theta[-1]), d[0], d[-1]], aspect='auto')
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f520d4347a0>
argsort(h.flat)[-3:]
array([32264, 32444, 27029])
for index in argsort(h.flat)[-3:]:
i_r, i_t = unravel_index(index,h.shape)
print(rad2deg(theta[i_t]), d[i_r])
-45.7541899441 -0.501385041551 -44.748603352 0.501385041551 -90.0 -29.5817174515
for index in argsort(h.flat)[-6:]:
i_r, i_t = unravel_index(index,h.shape)
print(rad2deg(theta[i_t]), d[i_r])
88.9944134078 32.5900277008 -88.9944134078 -28.5789473684 -88.9944134078 -27.5761772853 -45.7541899441 -0.501385041551 -44.748603352 0.501385041551 -90.0 -29.5817174515
# from http://stackoverflow.com/questions/9111711/get-coordinates-of-local-maxima-in-2d-array-above-certain-value
import scipy.ndimage.filters as filters
threshold = 0
neighborhood_size = 5
data_max = filters.maximum_filter(h, neighborhood_size)
maxima = (h == data_max)
data_min = filters.minimum_filter(h, neighborhood_size)
diff = ((data_max - data_min) > threshold)
maxima[diff == 0] = 0
imshow(maxima, aspect='auto', cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520d709fd0>
threshold = 30
neighborhood_size = 5
data_max = filters.maximum_filter(h, neighborhood_size)
maxima = (h == data_max)
data_min = filters.minimum_filter(h, neighborhood_size)
diff = ((data_max - data_min) > threshold)
maxima[diff == 0] = 0
imshow(maxima, aspect='auto', cmap=cm.gray, interpolation='nearest')
<matplotlib.image.AxesImage at 0x7f520d37f0d0>
argwhere(maxima)
array([[151, 0], [181, 45], [213, 178], [246, 90]])
for i_r, i_t in argwhere(maxima):
print(rad2deg(theta[i_t]), d[i_r])
-90.0 -29.5817174515 -44.748603352 0.501385041551 88.9944134078 32.5900277008 0.502793296089 65.6814404432
imshow(img, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520da52850>
from skimage import data
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.transform import warp, AffineTransform
from skimage.draw import ellipse
tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
translation=(210, 50))
image = warp(data.checkerboard(), tform.inverse, output_shape=(350, 350))
rr, cc = ellipse(310, 175, 10, 100)
image[rr, cc] = 1
image[180:230, 10:60] = 1
image[230:280, 60:110] = 1
imshow(image)
<matplotlib.image.AxesImage at 0x7f520bcb3810>
coords = corner_peaks(corner_harris(image), min_distance=5)
coords_subpix = corner_subpix(image, coords, window_size=13)
imshow(image, interpolation='nearest', cmap=plt.cm.gray)
plot(coords[:, 1], coords[:, 0], '.b', markersize=3)
plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15)
axis((0, 350, 350, 0));
coords = corner_peaks(corner_harris(lenabw), min_distance=5)
coords_subpix = corner_subpix(lenabw, coords, window_size=30)
imshow(lenabw, interpolation='nearest', cmap=plt.cm.gray)
plot(coords[:, 1], coords[:, 0], '.b', markersize=6)
plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=10)
axis((0, 512, 512, 0));
from skimage.filter import canny
coords = corner_peaks(corner_harris(canny(lenabw, 2, 0.05, 0.4)), min_distance=5)
coords_subpix = corner_subpix(lenabw, coords, window_size=30)
imshow(canny(lenabw, 2, 0.05, 0.4), interpolation='nearest', cmap=plt.cm.gray)
plot(coords[:, 1], coords[:, 0], '.b', markersize=6)
plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=10)
axis((0, 512, 512, 0));
/usr/lib/pymodules/python2.7/matplotlib/lines.py:503: RuntimeWarning: invalid value encountered in greater_equal return np.alltrue(x[1:] - x[0:-1] >= 0)
imshow(coins, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520bb84fd0>
coins_hist = measurements.histogram(coins, 0, 1, 256)
plot(linspace(0,1,256), coins_hist)
grid()
coins_hist = measurements.histogram(coins, 0, 1, 256)
plot(linspace(0,1,256), coins_hist)
grid()
vlines((0.39, 0.53, 0.6), 0, 1400)
<matplotlib.collections.LineCollection at 0x7f520b9dcd50>
imshow(coins > 0.39, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520bc0c710>
imshow(coins > 0.53, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520b91db50>
imshow(coins > 0.6, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520b85c690>
from scipy.ndimage.morphology import binary_dilation, binary_erosion
th = coins > 0.53
dilated = binary_dilation(th, ones((2,2)) )
imshow(dilated, cmap= cm.gray)
<matplotlib.image.AxesImage at 0x7f520b79c450>
eroded = binary_erosion(dilated, ones((4,4)) )
imshow(eroded, cmap= cm.gray)
<matplotlib.image.AxesImage at 0x7f520b6cff10>
th = coins > 0.5
dilated = binary_dilation(th, ones((2,2)) )
eroded = binary_erosion(dilated, ones((4,4)) )
imshow(eroded, cmap= cm.gray)
<matplotlib.image.AxesImage at 0x7f520b610c50>
th = coins > 0.5
from scipy.ndimage.morphology import binary_dilation, binary_erosion, binary_fill_holes
filled = binary_fill_holes(th)
imshow(filled, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520b550690>
eroded = binary_erosion(filled, ones((2,2)))
imshow(filled, cmap= cm.gray)
<matplotlib.image.AxesImage at 0x7f520b4920d0>
imshow(filled, cmap= cm.gray)
contour(filled, [0.5], linewidths=2, colors='r')
<matplotlib.contour.QuadContourSet instance at 0x7f520b4a03f8>
imshow(lenabw, cmap= cm.gray)
contour(lenabw, [0.7], linewidths=2, colors='r')
<matplotlib.contour.QuadContourSet instance at 0x7f520b293248>
from scipy import ndimage
label_objects, nb_labels = ndimage.label(filled)
sizes = bincount(label_objects.ravel())
imshow(coins, cmap= cm.gray)
contour(coins, [0.5], linewidths=2, colors='r')
<matplotlib.contour.QuadContourSet instance at 0x7f520ac7b248>
imshow(coins, cmap= cm.gray)
contour(contour, [0.5], linewidths=2, colors='r')
<matplotlib.contour.QuadContourSet instance at 0x7f520b165680>
print(sizes)
[78092 300 879 1 8 1 3 2 6 1 1 3 7 1 1 1 3 1 2 2 1 1 5 4 2 1 1 1 2 1 2 2587 3 2 1 1 4 1652 1598 1378 1160 1113 1 1 1 2 3 3 1 1833 1308 1198 1161 1107 1087 3 2949 1693 1457 1461 1044 1143 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2293 1348 1812 1697 1 1 1333 1 1419 1 3 3 24 1 1 2 18 1 2 2 1 1 1 11 1 3 1 1 1 1 1 2 1 1 1 2 3 3 1 1 4 3 5 2 4 2 1 1 1 5 1 1 3 1 1 2 1 4 1]
Threshold by size
mask_sizes = sizes > 1000
mask_sizes[0] = 0
coins_cleaned = mask_sizes[label_objects]
labeled_coins, nb_coins = ndimage.label(coins_cleaned)
imshow(labeled_coins)
<matplotlib.image.AxesImage at 0x7f520aaa0b90>
coins
array([[ 0.18431373, 0.48235294, 0.52156866, ..., 0.05490196, 0.01176471, 0.04705882], [ 0.36470589, 0.56470591, 0.56862748, ..., 0.04705882, 0.02745098, 0.02745098], [ 0.49411765, 0.57647061, 0.56078434, ..., 0.00784314, 0.05098039, 0.01176471], ..., [ 0.31764707, 0.30980393, 0.29019609, ..., 0.02352941, 0.01568628, 0.02745098], [ 0.34509805, 0.32156864, 0.29019609, ..., 0.01960784, 0.02745098, 0.03137255], [ 0.35686275, 0.30980393, 0.26666668, ..., 0.01568628, 0.03921569, 0.02745098]], dtype=float32)
from scipy.ndimage.measurements import center_of_mass
cmass = center_of_mass(coins_cleaned.T, labeled_coins.T, arange(nb_coins + 1))
/usr/lib/python2.7/dist-packages/scipy/ndimage/measurements.py:1277: RuntimeWarning: invalid value encountered in true_divide for dir in range(input.ndim)]
imshow(coins, cmap=cm.gray)
plot(*zip(*cmass), marker='o', ls= '')
[<matplotlib.lines.Line2D at 0x7f520acaba50>]
edges = gaussian_laplace(coins, 10) < 0
imshow(gaussian_laplace(coins, 10))
<matplotlib.image.AxesImage at 0x7f520a576590>
imshow(edges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520a937bd0>
eroded = binary_erosion(edges, ones((4,4)))
imshow(eroded, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520a881390>
eroded = binary_erosion(edges, ones((20,20)))
imshow(eroded, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520a7b5e90>
from scipy import ndimage
label_objects, nb_labels = ndimage.label(eroded)
sizes = bincount(label_objects.ravel())
imshow(label_objects)
<matplotlib.image.AxesImage at 0x7f520a6f4a10>
print(nb_labels)
24
from scipy.ndimage.measurements import center_of_mass
cmass = center_of_mass(eroded.T, label_objects.T, arange(nb_labels + 1))
imshow(coins, cmap=cm.gray)
plot(*zip(*cmass), marker='o', ls= '')
[<matplotlib.lines.Line2D at 0x7f520a69a4d0>]
from skimage import data, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
# Load picture and detect edges
edges = filters.canny(coins*255, sigma=3, low_threshold=10, high_threshold=50)
imshow(edges, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x7f520a4b7110>
# Detect two radii
hough_radii = np.arange(15, 30, 2)
hough_res = hough_circle(edges, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
num_peaks = 2
peaks = peak_local_max(h, num_peaks=num_peaks)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
centers
[array([201, 279]), array([ 58, 105]), array([56, 98]), array([127, 156]), array([198, 154]), array([124, 336]), array([266, 114]), array([125, 45]), array([193, 212]), array([ 51, 155]), array([118, 271]), array([262, 301]), array([260, 45]), array([261, 299]), array([ 44, 335]), array([258, 172])]
len(centers)
16
imshow(coins.T, cmap=cm.gray)
plot(*zip(*centers), marker='o', ls= '')
[<matplotlib.lines.Line2D at 0x7f520a47ced0>]
# Draw the most prominent 5 circles
image = coins[:]*255
image = color.gray2rgb(image)
for idx in np.argsort(accums)[::-1][:]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
image[cy, cx] = (220, 20, 20)
cx, cy = circle_perimeter(center_y, center_x, int(radius*1.1))
image[cy, cx] = (220, 20, 20)
imshow(image.astype(uint8), cmap=plt.cm.gray)
<matplotlib.image.AxesImage at 0x7f520a341f90>
fruits = imread('vitamy fruits.png')
imshow(fruits)
<matplotlib.image.AxesImage at 0x7f520a22ca50>
from skimage.segmentation import felzenszwalb, slic, quickshift
from skimage.segmentation import mark_boundaries
from skimage.util import img_as_float
img = img_as_float(fruits[:, :])
segments_fz = felzenszwalb(img, scale=100, sigma=1, min_size=50)
segments_slic = slic(img, n_segments=250, compactness=10, sigma=1)
segments_quick = quickshift(img, kernel_size=6, max_dist=6, ratio=0.5)
print("Felzenszwalb's number of segments: %d" % len(np.unique(segments_fz)))
print("Slic number of segments: %d" % len(np.unique(segments_slic)))
print("Quickshift number of segments: %d" % len(np.unique(segments_quick)))
Felzenszwalb's number of segments: 2693 Slic number of segments: 247 Quickshift number of segments: 2169
imshow(mark_boundaries(img, segments_fz))
title("Felzenszwalbs's method")
<matplotlib.text.Text at 0x7f5209b2f510>
imshow(mark_boundaries(img, segments_slic))
title("SLIC")
<matplotlib.text.Text at 0x7f5209ad3d90>
imshow(mark_boundaries(img, segments_quick))
title("Quickshift")
<matplotlib.text.Text at 0x7f5209a0bb50>
from scipy import ndimage
from skimage.morphology import watershed
from skimage.feature import peak_local_max
# Generate an initial image with two overlapping circles
x, y = np.indices((80, 80))
x1, y1, x2, y2 = 28, 28, 44, 52
r1, r2 = 16, 20
mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
image = logical_or(mask_circle1, mask_circle2)
imshow(image)
<matplotlib.image.AxesImage at 0x7f52098f1b90>
# Now we want to separate the two objects in image
# Generate the markers as local maxima of the distance to the background
distance = ndimage.distance_transform_edt(image)
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)),
labels=image)
markers = ndimage.label(local_maxi)[0]
labels = watershed(-distance, markers, mask=image)
subplot(131)
imshow(image, cmap=plt.cm.gray, interpolation='nearest')
title('Overlapping objects')
subplot(132)
imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
title('Distances')
subplot(133)
imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
title('Separated objects')
<matplotlib.text.Text at 0x7f5209777110>
By: Andrés Cabrera mantaraya36@gmail.com For MAT course MAT 201A at UCSB
This ipython notebook is licensed under the CC-BY-NC-SA license: http://creativecommons.org/licenses/by-nc-sa/4.0/
A lot of material adapted from the scikits image tutorials. http://scikit-image.org