Prune iteratively the leaves of a toy graph with graph-tool

A measure of the similarity of two weighted multigraphs?

To compare shapes, typically obtained after image segmentation, it has been shown that their skeleton could be used. A search with the keywords : skeleton+pruning+image+recognition yields a lot of results.

Assuming that the size of the branches of the skeleton, linking the end-points and the branched points together, can be represented by the weights of the edges in a graph, a weighted graph can represents naturally the skeleton of a shape.

From the skeleton of the B character (left) to its associated (networkx) graph (right).

Previously, pruning of skeletons was performed by removing pixels, one after an other, from the end-points by morphological hit-or-miss. A response curve of a skeleton to pruning, can be defined by counting the end-points of the skeleton as pruning is iterated. It was observed that the responses from similar shapes were similar. So, such response curves may be used to compute some similarity measures, or distances, between shapes.

The time necessary to prune a skeleton depends on its number of pixels. So it may be worth  prior pruning to convert a skeleton into a graph. Thus pruning a graph could be independent from the weight of its edges.

The response of a a weighted multi-graph to iterative pruning was performed with graph-tool (and additionnal python librairies : mahotas, numpy, scipy, matplotlib). This response yields a curve: a "Response to Pruning Curve" or RPC. Then a distance was defined between two RPCs, to compute a similarity index between two graphs.

The graph tool version used is:

2.11 (commit b6b2152c, Mon Oct 26 10:31:30 2015 +0100)

This post summarises the ideas developed in a jupyter notebook.

Algorithm for pruning a graph from its leaves :

1. Make a weighted un-directional multi-graph
2. (Visualize it)
3. Start pruning:
1. Count the vertices on degree 1
2. Read the weight of the edges bound to vertices of degree1 (V1)
3. Get the edge(s) of minimal weight
4. Remove this/these edge(s)
5. Decrease the weight of the other edges bound to vertices v1 by this minimal weight
4.  Do it until no more vertices of degree 1 remains
5. Cumulate the weight of removed edges
6. The response to pruning of a graph, associates the number of degree1 vertices to the cumulative weight of the removed edges.

Make a toy-graph to play with:

By default in graph-tool, the graphs are multi-graphs and that's what we need. The edges of the graph must be weighted, this property has to be created. The weight of an edge corresponds to the size of a branch in a skeleton (the number of pixels). The following python code, shows the creation of a graph with an edge property called 'weight': 

 

The graph can be displayed such that the thickness of its edges maps to the weight of edges:

Start to prune the graph:

Two different versions of the pruning function were written. The second called prune_once_faster() depends on the GraphView() method of the graph-tool library to filter vertices of degree 1, hoping to make a faster function:

The pruning function is destructive, applying it iteratively on the same graph yields, if you write in a jupyter cell:

G0 = make_toy_graph()
g0 = G0.copy()
prune_graph_once(g0)
g1 = g0.copy()
prune_graph_once(g0)
g2 = g0.copy()
prune_graph_once(g0)
g3 = g0.copy()

print G0, number_of_degree1_vertices(G0)
print g1, number_of_degree1_vertices(g1)
print g2, number_of_degree1_vertices(g2)
print g3, number_of_degree1_vertices(g3)

 You will have as output:

<Graph object, undirected, with 4 vertices and 4 edges at 0x9c65c1ac> 2
<Graph object, undirected, with 4 vertices and 3 edges at 0x9c65cc2c> 1
<Graph object, undirected, with 4 vertices and 2 edges at 0x9c654dec> 1
<Graph object, undirected, with 4 vertices and 1 edge at 0x9c65424c> 0

 Indicating that the number of edges decreases as the number of vertices of degree 1, each time the graph is pruned.

Pruning iteratively a graph: 

Let's count the number of vertices of degree 1 as the graph is pruned. The toy graph used here to check its response, has two vertices of degree 3 and two vertices of degree 1 (vertices 2 an 3). The initial values of the weight of the edges bound to vertices 2 and 3, are :

15, 20

The smallest weight is 15 So after one pruning, the weights becomes:

15-15 = 0, 20-15 = 5

The edge 1-2 of weight 0 is removed (vertex 2 is now of degree 0). 

The edge 1-3 has the smallest weight (5). 

So after a second round of pruning, the edge 1-3 is removed.

The vertex 1 was of degree 2, now it is of 1. The weight of the edge 0-1 is 8.

A third pruning removes the last edge 0-1 of weight equal to 8.

The degree of vertex 0 is 2 (there's a loop), so the edge 0-0 can't be pruned. The different weights are cumulated as follow:

0, 15, 15+5, 20+8

Applied to a skeleton made of pixels, this would mean that:

• 15 prunings of 1 pixel would remove one branche of size 15.

• 20 prunings would remove the two branches of size 15 and 20.

• 28 prunings would remove the three branches of sizes 15, 20 and 8 (initially between two vertices of degree 3).

The weights are stored in a vector X .

The corresponding counts of degree 1 vertices are stored in a vector Y.

Both X,Y are returned by the function response_to_iterative_pruning()

Make a feature vector from the pruning curve?

A response curve was plotted as a step curve:

The upper curve is the response when the weight of edges are used. In the lower one, the weights were normalised by the total weights of the graph (% total length or total weight), hopping to produce a scale invariant response. For the following, only the raw response will be used.

Distance between two graphs:

First, we need new toy graphs to compare them. The second new graph differe from the first only by its weight:

The third one has no loop:

 A fourth has a loop elsewhere:

The response of the three first graphs are:

Overlayed responses of first and third graph (right)

We could try to compute the area between two step curves. For two identical curves the area would be null, as the curves differe the area may increase. Fine!

However, the curves are just defined from discrete data, it may not be so easy to define for such kind of integration.

Let's forget the curves and just look at the points from the scatter plot. When looking at two pruning curves obtained from toy-graphs 2 and 3, we could just compute the distances matrix between the points from the two curves. The more the curves are different, the more the distances should increase ... why not, give it a try:

scipy provides just what we need!

The response curve is given as (X, Y) with

• X : The pruning values 

• Y : the corresponding number of leaf in the graph (vertices of degree 1)

First, let's compute the distances between the points of the same curve:

T=make_toy_graph()
X, Y = response_to_iterative_pruning(T.copy())
pruning_curve = zip(X,Y)
D = cdist( pruning_curve, pruning_curve)
print D

Thanks to scipy, we get a 4x4 (since the graph has four vertices) matrix of distances D :

[[  0.          15.03329638  20.02498439  28.0713377 ]
 [ 15.03329638   0.           5.          13.03840481]
 [ 20.02498439   5.           0.           8.06225775]
 [ 28.0713377   13.03840481   8.06225775   0.        ]]

Which can be displayed as an image (the matrix is symetric, only upper element are displayed):

For each of the three first graphs T,T1,T2 , one can define 32−32+3=6   matrices of distances.

Those matrices are elements of another matrix (a matrix of matrices). The elements of the diagonal are squared matrices corresponding to "intra-distances" (distances between the points of the same response curve).

The elements outside the diagonal correspond to "inter-distances" matrices (distances computed from points belonging to two different response curves), the "inter-distance matrices", here, can be  rectangular if the corresponding pairs of graphs do not have the same number of vertices (and square if the two graphs have the same number of vertices):

 For each matrix, the total distance can be computed as the sum of their elements. This yields:

[[ 178.  187.  256.]
 [   0.  105.  135.]
 [   0.    0.  140.]]

At first sight, this is not good, elements outside the diagonal (total distance between two response curves, that is the distance between two graphs)  should be greater than those of the diagonal.

Defining a normalized distance between two graphs

We try to heal this situation by normalising the inter-distances . The total inter-distances is normalized by the sum of total intra-distances. The total inter distances is multiplied by two in order to have a normalised distance:
 nD .  

nD(responsecurvei,responsecurvei)=1

When computed against the pair of same response curve, the normalized distance equals 1 as wished.

Finally nD  is computed as follow:

nD(ri,rj)=12distance(ri,rj)ri+rj=ri+rj2distance(ri,rj)

such that:

nD(ri,rj)<=1

Let's see the matrix of normalised distance on the three first graphs:

[[ 1.          0.75751403  0.62186779]
 [ 0.          1.          0.90636131]
 [ 0.          0.          1.        ]]

As an image:

This seems better. With that distance, the similarity of a graph on itself is 1 and two different graphs have a similarity bellow 1. It remains to see if that distance can discriminate between graphs of interest possibly after supervised classification.

An endpoints spectrum:possible features from the skeleton?

Looking for topological features:

The skeleton of a binary shape can be used to find features to classify it. The number of endpoints or the number of branched points should be good candidates. A shape which skeleton has no endpoint nor no branched point could be a circle:

 O

A shape with one branched point and no endpoint, may looks like this :

8

Binary model shapes were loaded in an ipython notebook with skimage which was capable of handling .pgm format :

Skeletons were generated by morphological thinning implemented in  mahotas:

endpoints (blue), branched points (green)

Endpoints Spectrum :

The skeleton of the hand has two endpoints close to top of each fingers, as it would be fine to have one endpoint per finger, let's prune the skeleton to meet this condition. This can be achieved after some trial, however if robust features could be found from the skeleton a systematic method is needed. So, let's prune the skeleton until it has no endpoint and count the number of pruning step necessary:

From the curve, the hand skeleton has 13 endpoints before any pruning, after 40 pruning steps, it has 6 endpoints. Around 90 pruning steps are necessary to remove all endpoints. By only looking at their curve, the hand shaped silhouette and the fish shaped silhouette are different.

Endpoints spectrum resists to scale variation after normalization:

As scale invariance is a nice feature, let's see what is happening when the scale of the shape is successively reduced by a factor of 0.9², 0.7², 0.5², 0.3²:

From the images above, it is clear that  the skeleton resists to scale variations but as the skeleton becomes smaller  less pruning steps are necessary to "burn" it to the last pruned pixels set. So let's norme the pruning steps by the initial skeleton size such for 0 pruning step, the normed pruning step is 100. As the size of the hand shaped silhouette decreases (from 81%, 49%, 25%, 9% of the original size), the endpoints spectrums  remain comparable:

Endpoints spectrum of similar shapes are similar:

How does the curve behave with similar shapes? Let's compare the endpoints spectrum of the serie of six hands:

The spectrums seem to resist to shape rotation (hand, hand-10°, hand90°) and to moderate shape deformation (hand-bent1,2). The sketch-hand silhouette is a right hand and its spectrum is comparable to the other spectrums which corresponds to left hands, so the endpoints spectrum is symetry invariant too.

Deriving a features vector from endpoints spectrum:

If shapes should be classified according to their topology, the whole endpoints spectrum is not necessary and only the transitions may be retained. Taking a  hypothetical skeleton whose endpoints spectrum before normalization is:

• For 7 prunings the skeleton has 0 endpoints

•       6                                        1

•       2                                        3

•       0                                        4

So a possible features vector would be:

(7,6,-1,2,0)

The value -1 indicates that the value 2 endpoints is missing. In the python script, the endpoints spectrum is coded as two lists:
             end_points_test   = [4,3,3,1,1,1,1,0]
             pruning_steps_test = [0,1,2,3,4,5,6,7]
A python function taking the steps list and the endpoints (spectrum) list yields a features vector :

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)
'''
     values = set(spectrum)
    feat_vect = [-1]*(len(values)+1)
    spec = np.asarray(spectrum)
    for epn in values:
        indexmin = np.where(spec == epn)[0][-1]
        feat_vect[epn] = steps[indexmin]
    return feat_vect
Supposing now the the skeleton size is 25, let's normalise the pruning steps as follow: 

 normed_pruning_steps_test = [100*(1-i/25.0) for i in pruntest] 

The endpoints spectrum is:

Let's take a normalised endpoints spectrum:

• normedsteps = [100.0, 96.0, 92.0, 88.0, 84.0, 80.0, 76.0, 72.0]

• endpoints = [4, 3, 3, 1, 1, 1, 1, 0] 

And let's submit it the previous function :

EPSpectrumToFeatures(normedsteps, endpoints)

 The function yields the vector:

[72.0, 76.0, -1, 92.0, 100.0]

which means that the skeleton has:

• 0 endpoint for a normed number of  prunings of 72%
• 1  ---------------------------------------------------------------------------76%
• no 2 endpoints, set normed number of  prunings to -1
• 3 -------------------------------------------------------------------- 92%
• 4 ------------------------------------------------------------------ 100%

The different endpoints spectrum curves suggest that the vector could used for classifying 2D binary shapes.

A feature vector from a real shape:

Let's consider the shape of a hand:

Its endpoints spectrum (normalised by skeleton size or not) was calculated and a feature vector was derived and displayed (red) with the spectrum:

Left: raw endpoints spectrum (blue) and its feature vector (red). Right: Normed endpoints spectrum and its feature vector, for 100% of the skeleton size, the skeleton has 14 endpoints, for ~94% the skeleton has 7 endpoints and for 83% of the skeleton size the skeleton has no endpoint.

For a larger number of image a python class, called Shape was written. It is instantiated as follow: 

shape= Shape(image)

The feature vector can be computed by:

vector = shape.End_Points_Features(normed = True)

or by:

vector = shape.End_Points_Features(normed = False)

Feature vectors from 216 images :

216 images of small 2D binary images (silhouettes) [Sharvit et al.] were used to derive feature vectors. This set of images consists in 18 type of silhouettes with 12 instance of each. For example, there are 12 birds silhouettes:

Prior any kind of classification the 216 normed feature vectors were aligned for visual inspection:

Some silhouettes, as the Bone silhouettes seems to show some intra-class structure, while Brick and Camel silhouettes seems different.

ipython notebook (updated 24/05/2013):

All the  tasks were performed in a ipython notebook available here.

 

Classifying chromosomes with scikit-learn

A trained dataset was build from 5 metaphases corresponding to 205 classified chromosomes or nuclei falling in one of four categories: single chromosome, overlapping chromosomes, nuclei or dusts(artefact of image segmentation).  

Two features were chosen to classify the particles:

• The particle area (normalised by image size)
• The convex-hull area of the particle ( area/convex hull). The idea is that particles that are clustered chromosomes will have a smaller ratio than single chromosomes.

The trained dataset is used to recognise 49 particles segmented from a sixth metaphase and then the predictions from the classifier compared to hand classified particles reached...

                                                             59.18%

# -*- coding: utf-8 -*-
"""
Created on Fri Oct  7 13:08:57 2011

@author: Jean-Patrick Pommier
"""
import KarIO
import os
import pylab
import pandas
import numpy as np
from sklearn import svm
#make a configurator object
config=KarIO.ClassifConf()
#build the the name feature file
##list all the feature files
##répertoire courant : os.listdir(os.getcwd())
featurespath=os.path.join(os.getcwd(),"Results","features")
labelpath=os.path.join(os.getcwd(),"Results","labels")
#print featurespath
##open features csv files
featuresFilesList=os.listdir(featurespath)
labelsFilesList=os.listdir(labelpath)
##
metalist=[0,1,2,3,4]
def makeDataShape(meta):
    featfile=config.user+'-'+config.slide+'-'+config.metaphases_list[meta]+'-'+config.counterstain+'.csv'
    labelfile='shapeCateg-'+featfile
    fea=pandas.read_csv(os.path.join(featurespath,featfile),header=None,index_col=None,names=['particle','ratio','area'])
    lab=pandas.read_csv(os.path.join(labelpath,labelfile),header=None,index_col=None,names=['particle','type'])
    del lab['particle']    
    #merge columns:features and label
    data=fea.join(lab)
    data.insert(0,'meta',int(config.metaphases_list[meta]))
    #print data
    return data

bigdata=makeDataShape(0)
for meta in metalist:
    bigdata=bigdata.append(makeDataShape(meta),ignore_index=True)

single=bigdata[bigdata['type']=='single']
touching=bigdata[bigdata['type']=='touching']
nuclei=bigdata[bigdata['type']=='nuclei']
dusts=bigdata[bigdata['type']=='dusts']

##ploting with pylab
#different colors according to the category

fig=pylab.figure()
ax = fig.add_subplot(111)
ax.scatter(single['ratio'],single['area'],c='green',marker='o')
ax.scatter(touching['ratio'],touching['area'],c='red',marker='o')
ax.scatter(nuclei['ratio'],nuclei['area'],c='blue',marker='o')
ax.scatter(dusts['ratio'],dusts['area'],c='pink',marker='o')

#train a classifier
trainedData=bigdata[bigdata['meta']<15]
untrained=bigdata[bigdata['meta']>=15]
print 'trained data'
print trainedData[:5]
#extract two columns from trainedData
#convert to numpy array
features=trainedData.ix[:,['ratio','area']].as_matrix(['ratio','area'])
test_features=untrained.ix[:,['ratio','area']].as_matrix(['ratio','area'])
print 'features'
print features[:5]
print 'features shape',features.shape
print 'features type',type(features)
##label is a string:single, touching,nuclei,dust
print 'labels convertion'
lab1=trainedData['type']
print 'lab1',type(lab1)
f=pandas.Factor(lab1)
print 'factor f',type(f)
print 'labels',f.labels[:5]
print 'labels type',type(f.labels)
print 'labels shape',f.labels.shape
#
##Classify with sklearn
classifier = svm.SVC()
model = classifier.fit(features,f.labels)
predicted=classifier.predict(test_features)

#match predicted /classified
hiddenlab1=untrained['type']
hiddf=pandas.Factor(hiddenlab1)
match=(predicted==hiddf.labels)
print"prediction"
print predicted[:5]
print 'true classification'
print hiddf.labels[:5]
print 'match'
print match[:5]
##Count sucess
success=np.sum(match[:]==True)
rate=100.0*success/(1.0*len(match))
print 'rate of good classification',success,'out of',len(match),'particles'
print rate,'% success'

chromosomes, overlapping chromosomes and nuclei

In an attempt to find features to distinguish chromosomes from overlapping chromosomes or from nuclei or small image segmentation artifacts ("dusts"), the area and the convex-hull area of isolated particles were computed after segmentation of a metaphasic chromosome image.

Image segmentation and convex-hull area computation were performed with scripts written previously. The features were saved in csv files.

miniKar produced an other csv file where the particle category (label) was saved.

This is an opportunity to handle the data with the pandas library and to see from this limited data set if the two features may be good candidates to classify at least particles in one of the three categories (single, overlapping, nuclei). A features file and the labels file were loaded, merged in one data frame. The features were filtered according to the label (the particle category) and then were displayed in a scatter plot:

green:single chromosomes; red overlapping chromosomes; blue:nuclei

The normalized area is the 1000*area of a particle divided by the image size (1024x1536). From this scatter plot, it's clear that at least an other feature will be necessary to separate the single chromosomes from the overlapping chromosomes. This may be due to the chromosomes bending.

A surprising thing, the ratio should be such: area/convexhull <1 and it is not true for all the particles.

import KarIO
import os,csv
import pylab
import pandas
#make a configurator object
config=KarIO.ClassifConf()
#build the the name feature file
##list all the feature files
##répertoire courant : os.listdir(os.getcwd())
featurespath=os.path.join(os.getcwd(),"Results","features")
labelpath=os.path.join(os.getcwd(),"Results","labels")
#print featurespath
##open features csv files
featuresFilesList=os.listdir(featurespath)
labelsFilesList=os.listdir(labelpath)
##
#let's find a pair of feature file/label file
#open with pandas
#two files match if they come from the same metaphase

featfile=config.user+'-'+config.slide+'-'+config.metaphases_list[2]+'-'+config.counterstain+'.csv'
labelfile='shapeCateg-'+featfile

fea=pandas.read_csv(os.path.join(featurespath,featfile),header=None,names=['particle','ratio','area'])
lab=pandas.read_csv(os.path.join(labelpath,labelfile),header=None,names=['particle','type'])
#merge columns:features and label
data=fea.join(lab)
single=data[data['type']=='single']
touching=data[data['type']=='touching']
nuclei=data[data['type']=='nuclei']
dusts=data[data['type']=='dusts']

##ploting with pylab
#different colors according to the category

fig=pylab.figure()
ax = fig.add_subplot(111)
ax.scatter(single['ratio'],single['area'],c='green',marker='o')
ax.scatter(touching['ratio'],touching['area'],c='red',marker='o')
ax.scatter(nuclei['ratio'],nuclei['area'],c='blue',marker='o')
pylab.show()