Shape analysis of 2D binary images by counting endpoints after successive skeleton pruning¶
image dataset from http://www.lems.brown.edu/vision/researchAreas/SIID/
Dependance: mahotas, pymporph, skimage, numpy¶
In [1]:
"""
Created on Sat Apr 20 12:57:37 2013
@author: Jean-Patrick Pommier
http://www.dip4fish.blogspots.fr
"""
#supposed to use mahotas 0.97
import mahotas as mh
from mahotas import polygon
import pymorph as pm
import numpy as np
from scipy import ndimage as nd
import skimage.transform as transform
import skimage.io as sio
import scipy.misc as sm
from matplotlib import pyplot as plt
import os
import math
In [2]:
def branchedPoints(skel):
branch1=np.array([[2, 1, 2], [1, 1, 1], [2, 2, 2]])
branch2=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 1]])
branch3=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 2]])
branch4=np.array([[2, 1, 2], [1, 1, 2], [2, 1, 2]])
branch5=np.array([[1, 2, 2], [2, 1, 2], [1, 2, 1]])
branch6=np.array([[2, 2, 2], [1, 1, 1], [2, 1, 2]])
branch7=np.array([[2, 2, 1], [2, 1, 2], [1, 2, 1]])
branch8=np.array([[2, 1, 2], [2, 1, 1], [2, 1, 2]])
branch9=np.array([[1, 2, 1], [2, 1, 2], [2, 2, 1]])
br1=mh.morph.hitmiss(skel,branch1)
br2=mh.morph.hitmiss(skel,branch2)
br3=mh.morph.hitmiss(skel,branch3)
br4=mh.morph.hitmiss(skel,branch4)
br5=mh.morph.hitmiss(skel,branch5)
br6=mh.morph.hitmiss(skel,branch6)
br7=mh.morph.hitmiss(skel,branch7)
br8=mh.morph.hitmiss(skel,branch8)
br9=mh.morph.hitmiss(skel,branch9)
return br1+br2+br3+br4+br5+br6+br7+br8+br9
def endPoints(skel):
endpoint1=np.array([[0, 0, 0],[0, 1, 0],[2, 1, 2]])
endpoint2=np.array([[0, 0, 0],[0, 1, 2],[0, 2, 1]])
endpoint3=np.array([[0, 0, 2],[0, 1, 1],[0, 0, 2]])
endpoint4=np.array([[0, 2, 1],[0, 1, 2],[0, 0, 0]])
endpoint5=np.array([[2, 1, 2],[0, 1, 0],[0, 0, 0]])
endpoint6=np.array([[1, 2, 0],[2, 1, 0],[0, 0, 0]])
endpoint7=np.array([[2, 0, 0],[1, 1, 0],[2, 0, 0]])
endpoint8=np.array([[0, 0, 0],[2, 1, 0],[1, 2, 0]])
ep1=mh.morph.hitmiss(skel,endpoint1)
ep2=mh.morph.hitmiss(skel,endpoint2)
ep3=mh.morph.hitmiss(skel,endpoint3)
ep4=mh.morph.hitmiss(skel,endpoint4)
ep5=mh.morph.hitmiss(skel,endpoint5)
ep6=mh.morph.hitmiss(skel,endpoint6)
ep7=mh.morph.hitmiss(skel,endpoint7)
ep8=mh.morph.hitmiss(skel,endpoint8)
ep = ep1+ep2+ep3+ep4+ep5+ep6+ep7+ep8
return ep
def pruning(skeleton, size):
'''remove iteratively end points "size"
times from the skeleton
'''
for i in range(0, size):
endpoints = endPoints(skeleton)
endpoints = np.logical_not(endpoints)
skeleton = np.logical_and(skeleton,endpoints)
return skeleton
def CountEndPoints(skel):
ep = endPoints(skel)
lab, n = mh.label(ep)
return n
def EndPointsSpectrum(skel):
#if no end points: its a closed curved :end
#count end points
#prune skeleton
#do it again until skeleton vanish, or has no endpoints
spectrum = []
steps = []
iteration = 0
epn = CountEndPoints(skel)
skellab, n = mh.label(skel)
skelsize = mh.labeled.labeled_size(skellab)
#print 'skelsize',skelsize[1]
if epn == 0 :
spectrum.append(epn)
return steps, spectrum
while skelsize[1] >= 1:
skel = pruning(skel, 1)
epn = CountEndPoints(skel)
if epn == 0 : #loop(s) in the skeleton
steps.append(iteration)
spectrum.append(epn)
return steps, spectrum
spectrum.append(epn)
steps.append(iteration)
iteration = iteration + 1
skellab, n = mh.label(skel)
skelsize = mh.labeled.labeled_size(skellab)
#print iteration, skelsize[1], epn
return steps, spectrum
def normedEPS(binImg):
sk = mh.thin(binImg)
lab, n = mh.label(sk)
size = mh.labeled.labeled_size(lab)
#print size, n
size = size[1]
steps, spectrum = EndPointsSpectrum(sk)
normSteps = [100-100.0*s/size for s in steps]
return normSteps, spectrum
def CountBranchedPoints(bIm):
sk = mh.thin(bIm)
bp = branchedPoints(sk)
lab, n = mh.label(bp)
return n
def CountHoles(bIm):
noHoles = mh.close_holes(bIm)
anti = np.logical_not(bIm)
holes = np.logical_and(noHoles, anti)
lab, n = mh.label(holes)
return lab,n
def segmentHolesBellow(bimg, thresh=5):
bimg = bimg>0
#filled convexhull
cvxf = polygon.fill_convexhull(bimg)
area= np.sum(bimg)
cvxf_areaF = np.sum(cvxf > 0)
#print area_wheel, areaF_wheel, 100.0*area_wheel/areaF_wheel
#holes size
holes_map, nbh = CountHoles(bimg)
holes_size_map = np.copy(holes_map)
holes_size_map = np.float32(holes_size_map)
for hole in range(1, nbh+1):
hole_size = np.sum((holes_map==hole)>0)
holes_size_map[holes_size_map==hole]=100.0*hole_size / 1.0*areaF_wheel
holes_size_map[holes_size_map>=thresh] = 0
smallHoles = holes_size_map>0
smallHolesFilled = np.logical_or(wheel<150, smallHoles)
return smallHolesFilled
def EPSpectrumToFeatures(steps, spectrum):
'''
eptest = [4,3,3,1,1,1,1,0]
pruntest = [0,1,2,3,4,5,6,7]
the returned vec:[7,6,-1,2,0] which means:
0 enpoints -> 7 prunings necessary
1 endpoint -> 6 prunings
4 endpoints-> 0 pruning (initial skeleton size)
'''
#feat_vect = [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]
values = set(spectrum)
feat_vect = [-1]*(len(values)+1)
spec = np.asarray(spectrum)
print 'EPS spec to feat func'
print 'values',values
print spec
print feat_vect
for epn in values:
indexmin = np.where(spec == epn)[0][-1]
feat_vect[epn] = steps[indexmin]
#print epn,np.where(spec == epn),np.where(spec == epn)[0],np.where(spec == epn)[0][-1], feat_vec
return feat_vect
def step(endpoints):
'''ep = [8,7,4,4,3,3,3,2,0]
'''
values = set(endpoints)
feat = [-1]*(max(values)+1)
endpoints = np.asarray(endpoints)
#print endpoints
#print values
for ep in values:
elements = np.where(endpoints == ep)[0]
last = elements[-1]
feat[ep] = last
return feat
#essai eps to feat
eptest = [4,4,3,3,1,1,1,1,0]
pruntest = [0,1,2,3,4,5,6,7,8]
normedprun = [100*(1-i/25.0) for i in pruntest]
print 'normesSteps', normedprun
print 'feat vector',EPSpectrumToFeatures(normedprun,eptest)
##
### Essai2
p1 = [0, 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]
ep =[13, 13, 12, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,6, 6,6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2,2,2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 0]
print len(ep)
print len(p1)
#print 'feat vector',EPSpectrumToFeatures(p1, ep)
#print 'normed feat vector',EPSpectrumToFeatures(normedprun,eptest)
#structuring elements
disk7 = pm.sedisk(3)#size 7x7: 7=3+1+3
disk5 = pm.sedisk(2)
normesSteps [100.0, 96.0, 92.0, 88.0, 84.0, 80.0, 76.0, 72.0, 68.0] feat vector EPS spec to feat func values set([0, 1, 3, 4]) [4 4 3 3 1 1 1 1 0] [-1, -1, -1, -1, -1] [68.0, 72.0, -1, 88.0, 96.0] 88 88
particle is a class to analyse particle anisotropy¶
In [3]:
class particle(object):
def __init__(self,arrayIm):
self.particuleImage=arrayIm
def orientation(self,step):
""" Performs successive rotations of a binary particle from 0 to 180 by "step"
Then at each angle performs:
*successives erosions
*pixel counts
Find relationship between angle and the minimal erosion to destroy the particle
An anisotropic particle such a long chromosome should have a
principal orientation, overlapping chromosomes should have at least two
orientation.
"""
def erodeVParticle(im):
'''Counts the number of vertical erosion neccessary
to destroy a particule
vline is a 3X1 structuring element
DestroyParticle returns the scalar number:nb_erosion
necessary to erode totally a binary particle'''
def countTruePix(bim):
'''bim must be a flat binary (True,False) array
counts the occurence of True (pixel) in a binary array/image '''
return np.sum(bim[:,:]==True)
vline=np.array([[0,1,0],
[0,1,0],
[0,1,0]])
#threshold the image
binaryIm=(im>0)
pixcount=countTruePix(binaryIm)
nb_erosion=1
while pixcount>0:
erodedIm=nd.binary_erosion(binaryIm,structure=vline,iterations=nb_erosion)
nb_erosion=nb_erosion+1
pixcount=countTruePix(erodedIm)
#erodedParticle work should be done here
return nb_erosion
#
#Analysis of the particle orientation starts here
#
maxRotationNubr=180/step
compassTable=np.zeros((2,maxRotationNubr+1),dtype=np.uint8)#rotation angle from 0 to 180 by step of 5°
#Store
#print compassTable.shape
i=0
angle=i*step
rotatedIm=self.particuleImage
#print self.particuleImage.shape
while i<=maxRotationNubr:
maxErosion=erodeVParticle(rotatedIm)
compassTable[0,i]=angle
compassTable[1,i]=maxErosion
#print "i=",i," angle=",angle," -erosion max:",maxErosion
rotatedIm=nd.rotate(self.particuleImage,angle)
i=i+1
angle=i*step
return compassTable
Let's use a class to handle shapes and to produce a feature vector¶
In [15]:
class Shape(object):
def __init__(self, silhouette, neg = False):
if neg == True:
self.silhouette = np.logical_not(silhouette > 0)
if neg == False:
self.silhouette = silhouette > 0
self.skeleton = mh.thin(self.silhouette)
self.skelsize = np.sum(self.skeleton > 0)
self.endPointSpectrum = [EndPointsSpectrum(self.skeleton)]
self.prunedSkeletons = {}
self.prunedSkeletons[0] = self.skeleton
def End_Points(self):
return endPoints(self.skeleton)
def Branched_Points(self):
return branchedPoints(self.skeleton)
def End_Points_Spectrum(self):
'''From skeleton image return end points spectrum
'''
return EndPointsSpectrum(self.skeleton)
def Normed_EPS(self):
return normedEPS(self.skeleton)
def End_Points_Features(self, normed = True):
''''return transitions of the endpoints spectrum as the number of
pruning steps if normed = False or
as a % of skeleton size if normed = True (100*(1-Pn/skelsize), n=0,1,...).
-1 value(s) are missing endpoints value for a given bumber of pruning
'''
eps = self.End_Points_Spectrum()
#neps = h0.NormedEPS()
feat = step(eps[1])
if normed == False:
#print 'raw'
return feat
elif normed == True:
#print 'normed'
size = self.skelsize
nfeat =[100*(1-f/(1.0*size)) for f in feat]
nfeat =np.asarray(nfeat)
nfeat[np.where(nfeat>100)]=-1
return nfeat
def feat(self):
print type(self.endPointSpectrum)
endpoints = self.endPointSpectrum[1]
return step(self, endpoints)
Load images¶
In [5]:
hands = ['hand','hand-10','hand-90','hand-bent1','hand-bent2','sketch-hand']
user=os.path.expanduser("~")
workdir=os.path.join(user,"Images","silhouette-database","hands")
handsImg = []
for img in hands:
image = sio.imread(os.path.join(workdir,img+'.pgm'))
handsImg.append(np.logical_not(image>0))
i = 1
title = iter(hands)
print len(handsImg)
for im in handsImg:
plt.subplot(int('16'+str(i)),frameon=False, xticks=[], yticks=[])#
plt.title(title.next(), fontsize = 8)
plt.imshow(im>0)
i = i+1
#plt.imshow(handsImg[0])
6
In [6]:
fishes = "/home/simon/Images/silhouette-database/fish"
print os.listdir(fishes)
fishesImg=[]
for im in os.listdir(fishes):
image = sio.imread(os.path.join(fishes,im))
fishesImg.append(np.logical_not(image>0))
plt.imshow(fishesImg[2])
['sketch.pgm', 'puffers.pgm', 'cardinalfishes.pgm', 'dogfish-sharks.pgm', 'goatfishes.pgm', 'whale-sharks.pgm', 'livebearers2.pgm', 'mullets.pgm', 'bonefishes.pgm', 'swordfishes.pgm', 'herrings.pgm', 'sketch2b.pgm', 'sand-tigers.pgm']
Out[6]:
<matplotlib.image.AxesImage at 0xbb4f7ec>
Detecting end-points or branched-points after skeleton pruning¶
In [7]:
if __name__ == "__main__":
handsSkel = []
#
prunedHandskel = {}
fishesSkel = []
# Skeleton+Pruning
for im in handsImg:
skel = mh.thin(im)
handsSkel.append(skel)
for im in fishesImg:
skel = mh.thin(im)
fishesSkel.append(skel)
hand0 = Shape(handsImg[0])
#Example with a hand
skel = hand0.skeleton
endpts = hand0.End_Points()
endpts = mh.morph.dilate(endpts, disk5)
brpts = hand0.Branched_Points()#
brpts = mh.morph.dilate(brpts, disk5)
#Example with a fish
fish0 = Shape(fishesImg[2])
skel_f = fish0.skeleton
endpts_f = fish0.End_Points()
endpts_f = mh.morph.dilate(endpts_f, disk5)
brpts_f = fish0.Branched_Points()#
brpts_f = mh.morph.dilate(brpts_f, disk5)
display_H0 = pm.overlay(hand0.silhouette, red = skel>0, blue = endpts > 0, green = brpts > 0)
display_F0 = pm.overlay(fish0.silhouette, red = skel_f>0, blue = endpts_f > 0, green = brpts_f > 0)
plt.subplot(121)
plt.imshow(display_H0)
plt.subplot(122)
plt.imshow(display_F0)
Comparing the evolution end-points number of a skeleton with pruning steps in two shapes (endpoints spectrum)¶
In [8]:
plt.subplot(111)
S= hand0.endPointSpectrum
pruning_steps= S[0][0]
end_points = S[0][1]
F= fish0.endPointSpectrum
pruning_steps_f= F[0][0]
end_points_f = F[0][1]
#print S[0]
#plt.gca().invert_xaxis()
plt.step(pruning_steps, end_points)
plt.scatter(pruning_steps, end_points)
plt.step(pruning_steps_f, end_points_f)
plt.xlabel('pruning steps')
plt.ylabel('end-points')
plt.legend(['hand','fish'])
plt.show()
Effect of image size on endpoints spectrum stability¶
Reduce image size successively: 100%, 81%, 49%, 25%, 9% (from initial size). The curve is more or less stable as the image resolution decreases.¶
In [9]:
im100 = pm.asf(handsImg[0], seq='CO', n=0)
scales =[1, 0.9, 0.7, 0.5, 0.3]
images =[]
for s in scales:
images.append(transform.rescale(im100, (s,s)))
#images.append(sm.imresize(im100, (s,s), interp='bilinear', mode=None)
i = 1
title = iter(scales)
for im in images:
plt.subplot(int('15'+str(i)),frameon=False, xticks=[], yticks=[])
plt.title(title.next())
display = pm.overlay(im>0, red = mh.thin(im>0)>0)
plt.imshow(display)
i = i+1
print len(images)
5
Show the end-points spectrum¶
In [10]:
scales =[1, 0.9, 0.7, 0.5, 0.3]
for im in zip(images, scales):
plt.xlim(75, 101)
plt.gca().invert_xaxis()
normstp, spec = normedEPS(im[0]>0)
plt.title('End-points Spectrum')
plt.step(normstp, spec)
#print 'scale:',im[1],' feat vect:',EPSpectrumToFeatures(normstp, spec)
#a = '100\left(1-\frac{pruning}{skel \; size} \right)'
plt.xlabel(r'normed pruning steps $ 100 \left(1-\frac{pruning}{skel \; size} \right)$', fontsize=16)#
plt.ylabel('end-points')
plt.legend(scales, title="scale factor")
plt.show()
comparison of several hands endpoints spectrum¶
In [11]:
hands = ['hand','hand-10','hand-90','hand-bent1','hand-bent2','sketch-hand']
for h in handsImg:
hand = Shape(h)
plt.xlim(75, 101)
plt.gca().invert_xaxis()
spec = hand.Normed_EPS()
normstp = spec[0]
ep = spec[1]
plt.title('End-points Spectrum')
plt.step(normstp, ep)
plt.legend(hands)
plt.show()
Extraction of feature vector¶
In [16]:
plt.imshow(handsImg[0])
h0 = Shape(handsImg[0])
pr = h0.End_Points_Spectrum()[0]
ep = h0.End_Points_Spectrum()[1]
print h0.skelsize
eps = h0.End_Points_Spectrum()
neps = h0.Normed_EPS()
feat = step(eps[1])
nfeat =[100*(1-f/(1.0*h0.skelsize)) for f in feat]
nfeat =np.asarray(nfeat)
nfeat[np.where(nfeat>100)]=-1
plt.subplot(121)
plt.xlim((-2,100))
plt.step(eps[0],eps[1])
plt.scatter(eps[0],eps[1])
plt.scatter(feat, [-0.5]*len(feat),c='red',marker='o')
plt.subplot(122)
plt.xlim((80,101))
plt.gca().invert_xaxis()
plt.step(neps[0],neps[1])
plt.scatter(nfeat, [-0.5]*len(nfeat),c='red',marker='o')
plt.scatter(neps[0],neps[1])
490
Out[16]:
<matplotlib.collections.PathCollection at 0xb83a44c>
In [17]:
print h0.End_Points_Features(normed=False)
print h0.End_Points_Features(normed=True)
[87, 86, 85, 69, 59, 57, 54, 30, 12, 4, -1, -1, 2, 1] [ 82.24489796 82.44897959 82.65306122 85.91836735 87.95918367 88.36734694 88.97959184 93.87755102 97.55102041 99.18367347 -1. -1. 99.59183673 99.79591837]
Load silhouette images¶
In [27]:
from collections import defaultdict
import skimage.io as io
silhouette216 = ['bird','Bone','brick','camel','car','children','classic','elephant','face',
'fork','fountain','Glas','hammer','Heart','key','Misk','ray','turtle']
user=os.path.expanduser("~")
workdir=os.path.join(user,"Images","shapes216")
filesname =defaultdict(list)
for files in os.listdir(workdir):
for silhouette in silhouette216:
if files.startswith(silhouette):
filesname[silhouette].append(files)
#print filesname['fork']
images_collections =defaultdict(list)
for silhouette in silhouette216:
collection = io.ImageCollection(os.path.join(workdir,silhouette+'*'))
print silhouette,':%s images added'%len(collection),collection[0].shape, collection[0].dtype
images_collections[silhouette].append(collection)
print type(images_collections['bird'])
print type(images_collections['bird'][0]), len(images_collections['bird'][0])
print type(images_collections['bird'][0][0])
bird :12 images added (108, 115) uint8 Bone :12 images added (127, 141) uint8 brick :12 images added (57, 113) uint8 camel :12 images added (108, 115) uint8 car :12 images added (63, 119) uint8 children :12 images added (123, 126) uint8 classic :12 images added (52, 115) uint8 elephant :12 images added (115, 115) uint8 face :12 images added (130, 109) uint8 fork :12 images added (42, 115) uint8 fountain :12 images added (90, 109) uint8 Glas :12 images added (100, 100) uint8 hammer :12 images added (47, 115) uint8 Heart :12 images added (95, 113) uint8 key :12 images added (60, 116) uint8 Misk :12 images added (129, 129) uint8 ray :12 images added (100, 138) uint8 turtle :12 images added (80, 115) uint8 <type 'list'> <class 'skimage.io.collection.ImageCollection'> 12 <class 'skimage.io._io.Image'>
Make a feature vector for each silhouette¶
In [30]:
feature_vectors_collection =defaultdict(list)
for silhouette in silhouette216:
print silhouette
for image in images_collections[silhouette][0]:
shape = Shape(image)
feature_vectors_collection[silhouette].append(shape.End_Points_Features(normed = True))
print type(feature_vectors_collection['bird'])
print type(feature_vectors_collection['bird'][0]), len(feature_vectors_collection['bird'][0])
print type(feature_vectors_collection['bird'][0][0])
bird Bone brick camel car children classic elephant face fork fountain Glas hammer Heart key Misk ray turtle <type 'list'> <type 'numpy.ndarray'> 12 <type 'numpy.float64'>
In [60]:
#f = feature_vectors_collection['bird'][0]
#plt.scatter(f, [-0.5]*len(f),c='red',marker='o')
import itertools
y = 0
#silh = reversed(silhouette216)
colors = itertools.cycle(["r", "b", "g","c","m","y"])
plt.figure(num=None, figsize=(8, 16), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(111,yticks=[])
plt.ylim((-1,220))
plt.xlim((85,102))
plt.gca().invert_xaxis()
plt.title('Aligned feature vectors from 18 family of shapes')
h=1
for silhouette in silhouette216:
n=1
col=next(colors)
#plt.legend(silhouette216, loc='center left', bbox_to_anchor=(1, 0.5))
#print silhouette, 6+12*(h-1)
plt.text(101.5,6+12*(h-1), silhouette)
h = h+1
for feat in feature_vectors_collection[silhouette]:
plt.scatter(feat, [y]*len(feat),s=8,c=col,marker='x')
n=n+1
y = y+1
print y, y/12.0
216 18.0
In [57]:
birds = images_collections['bird'][0]
print type(birds), len(birds),type(birds[0])
i=1
#for l in range(1,4):
# for c in range (1,5):
# plt.subplot(l,c,i)
# plt.imshow(birds[i-1])
# i = i+1
index = [i for i in range(1,13)]
for bird, im in zip(birds, index):
plt.subplot(3,4,im, frameon=False, xticks=[], yticks=[])
plt.imshow(bird<128)
<class 'skimage.io.collection.ImageCollection'> 12 <class 'skimage.io._io.Image'>
In [ ]: