version 0.2 24/05/2013 """ 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 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) 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 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 shape images (pgm format) 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]) load fish images 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]) 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) 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() 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) 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() 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() 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]) print h0.End_Points_Features(normed=False) print h0.End_Points_Features(normed=True) 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]) 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]) #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 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)