Hunting blobs

 Pictures from my last safari:

Graph again

Update 06 /11 /2014

The ipython notebook from the previous post was rewritten . A python  function building a graph from a skeleton was modified to generate a multigraph (networkx). This function (available in an ipython notebook) is used as follow:

    Graph= nx.MultiGraph()
    C8_Skeleton_To_Graph_01(Graph, skeleton)

With:

• Graph : a networkx multigraph object.
• skeleton : the binary image of skeleton obtained by morphological thining with mahotas.

In the following example (from the previous ipython notebook), the skeleton of the letter 'a' was converted into a networkx multigraph: 

image_test= makeLetterImage('a', 75)
skeleton_test = mh.thin(image_test)
_,_,Bp_test,Ep= SkeletonDecomposition(skeleton_test)
Ep_Ep = skeleton_test*np.logical_not(Ep)
Ep_Ep,_ = mh.label(Ep,Ep)
l_Ep, _ = mh.label(Ep)


Graph_test = nx.MultiGraph()
C8_Skeleton_To_Graph_01(Graph_test, skeleton_test)
print Graph_test.edges(data=True)

figsize(16,8)
subplot(131)
imshow(skeleton_test+l_Ep+3*Bp_test,interpolation='nearest')
subplot(132, xticks=[],yticks=[])
nx.write_dot(Graph_test,'multi.dot')
!neato -T png multi.dot > multi.png
imshow(mh.imread('multi.png'), interpolation='nearest')
subplot(133)
nx.draw(Graph_test)

The skeleton has two edges, linking the same branched points (yellow and red pixels). The graph was exported to build an image of the multigraph (middle), networkx doesn't seem to be capable to dispolay correctly such multigraph (right)

When applied to a lowercase alphabet:

The function produced the corresponding graphs (the letter 'o' excepted):

However,it is possible to find case where the function failed to buil a graph as:

   imTest = makeLetterImage('b', 70)
   skelTest = mh.thin(imTest)
   Ep_Bp, Bp_Bp, Bp, Ep = SkeletonDecomposition(skelTest)

where:

• skeltest: the image skeleton obtained by thining with mahotas.
• Bp : image of branched-point(s)
• Ep: image of end-point(s)
• Bp_Bp :edge(s) between branched-points
• Ep_Bp:edge(s) between end-point(s) and branched-point(s)

If the skeleton of the letter 'b' is decomposed as follow:

The function fails to produce a graph due to the presence of a closed edge (right) which should be processed as a self loop for the corresponding branched-point (middle image).

An other graph library: Graph-tool

The graphic representation of a multigraph produced by graphtool is so nice that it may worth to learn the graphtool syntax to make a graph:

Construct a graph from the skeleton image of a binary form

Previously an attempt was made to write a code capable of building a graph from the skeleton of an image and it was not a complete success since some preprocessing (pruning) had to be done to build a graph. Here an other python code was written and at very first sight, it seems to produce a graph without complaining....sometimes (see at the end).

The code was run in a ipython notebook running on sagemath cloud. It relies on several libraries : mahotas and networkx, numpy is provided by default. Scikit image was also tested. Pillow was used to generate images of some letter as model.

The idea is the following:

• make the skeleton of a binary image.
• get its branched-points and add them in a graph.
• get its end-points, add them to the graph.
• get the edges of the skeleton and decompose them into:
• edges between end-points and branched-points.
• edges between branched points.
• add edges to the graph between end-points and branched-points.
• add edges to the graph between branched-points.

Images models

Since the code is run on sagemath cloud, Pillow has to be installed. From a terminal :

pip install Pillow

 The verdana true type font was obtained from fontsupply and the file verdana.ttf had to be at the same location than the ipython notebook ( verdana.ttf was uploaded on SMC project at the same location than the ipython notebook).

The code used to generate images and draw them is:

and we get:

What kind of skeleton?

Both mahotas an scikit-image implement morphological operators to generate the skeleton of a binary shape. Scikit-image seems to use a neighbourhood  on 4 pixels (C4) and mahotas on on neighbourhood of 8 pixels (C8). This is why mahotas was retained to produce skeletons:

• medial axis  of the letter A computed with scikit image (left), branched points and end-points are extracted using hit-or-miss operator implemented in mahotas (middle). The branched-points can be removed from the medial axis and the resulting image can be labelled (left), pixels connected only vertically or horizontally share the same label.

Skeleton from the letter A (left)

•  If a skeleton is generated by morphological thining with mahotas, we get the left image. The fourth image is the skeleton deprived of its branched points. Scikit-image provide a thining operator with a C4 neighborhood (right image):

 

C8 skeleton by thining (left). C4 skeleton by thining (right).

What's happen if a C8 skeleton is labeled? Here, both mahotas and scikit-image provide the same result:

edges (mahotas:left; scikit-image:right)

The labelled regions fit here what can be seen as an edges.

Prior to start to build a graph, several additionnal images are needed:

Edges decomposition:

Let's decompose the image of edges by considering edges connecting end-points to branched-points and edges connecting only branched points:

Starting to build the graph by the branched-points:

 The branched-points are added first. The branched point of label 1 in the image, is the first vertex in the graph, and so on:

The end points are the added to the graph, there is a shift between the end-points label in the image and their index in the graph. The first end-point in the image will be the fith vertex in the graph:

The vertex of the graph (G) have attributes derived from the images they come from: their position 'pos' define by (row,col), their label in the image of origin. Typing in a ipython cell:

          print G.nodes(data=True)
          print nx.get_node_attributes(G,'label')

yields:

[(1, {'pos': (3, 9), 'label': 1}), (2, {'pos': (3, 12), 'label': 2}), (3, {'pos': (17, 5), 'label': 3}), (4, {'pos': (17, 17), 'label': 4}), (5, {'pos': (1, 9), 'label': 1}), (6, {'pos': (1, 13), 'label': 2}), (7, {'pos': (25, 3), 'label': 6}), (8, {'pos': (25, 21), 'label': 8})]
{1: 1, 2: 2, 3: 3, 4: 4, 5: 1, 6: 2, 7: 6, 8: 8}

Linking end-points to branched points:

 The idea here is to get the neighbourhood of each branched-point in the "edges ep-ep" image and to add an edge in the graph between the corresponding branch-point and end-point(s). This yield:

Linking branched-points

•  A dictionnary where the keys are the labels of the branched-points and the values are the label of the "edges bp-bp" in the neighbourhood of the branched point is build. Thus there's one key (bp label) for several values (edges labels).
• But an edge can be seen as one key (the edge label) for two values (the two labels of the two branched-points). So the dictionnary has to be inverted in that very particular way. Fortunately, that problem was solved on stackoverflow.

 Finaly, we have:

 the whole steps are gathered in a function:

Processing other  letters gives for example:

The ipython notebook is here:

View the notebook on ipython.org

Sometimes building the graph fails:

On a larger set of images:

• For the image corresponding to the letter K, the graph construction fails due to the presence of branched regions made of several pixels instead of one.
• There's also a problem with the letter D, its graph has no loop.
• Since there's no branched point in the skeleton of the letter C, its graph is not build.

 The ipython notebook has been reordered and is available on ipython.org

Resolving overlapping chromosomes: an emblematic case

Metaphasic chromosomes preparations contain most of the time overlapping chromosomes as follow :

A cluster of two overlapping chromosomes

Here a cluster of two chromosomes was isolated and the two chromosomes were segmented using a contour based method implemented in an ipython notebook.

First a polygonal approximation of the contour is performed and negative corners are extracted as previously:

Polygonal approximation of the cluster contour (red line). Negative corners (yellow dots), positive corners (red dots).

Four points (the yellow ones) are then used to define a quadrilateral of maximal area:

The quadrilateral is converted to an image (overlapp domain:right image) and substracted from the original thresholded image yielding here four chromosomal elements (left image):

The four chromosomal elements are labelled:

Left:labelled chromosomal elements (the size is figured above the image). Right: labelled overlapping domain.

and isolated:

The idea is to used combinations of elements to build candidate chromosomes. Here the image is supposed to contain two chromosomes, so combinations of two elements in a set of four elements yield the following candidates:

For each candidate, several descriptors are computed: the number of connected components (cc), the eccentricity (ecc), the ratio chromosome area over its convex hull area (r) and the number of negative corner on the contour.

The two first highest ratio correspond to the two overlapping chromosomes which can be isolated:

Download the ipython notebook from google drive

Or view it

Which structuring elements to detect branched points?

In a 2D skeleton obtained after thinning of a binary image, pixels having 3 or 4 (or more...) neighbors are branched-points. Hit-or-miss operator can be used to detect branched points (or junction points). Structuring elements, SE, can be 3x3 arrays representing the pixels configuration to be detected. Here different 3x3 structuring elements were explored on some simple test images with the hit-or-miss operator implemented in mahotas 1.0 .

• There are two SE to detect pixels with 4 neighbors, called here X0 and X1.
• There are two set of SE detecting pixels with at least 3 pixels, a Y-like configuration and a T-like configuration. Each of these configurations can be rotated by 45°, making 16 different SE.

18 SE were found , since don't care pixels were used a branched points can be detected by different SE, a pixel detected by X0 should be detected by T4 too. The 18 SE can be displayed as 3x3 image:

Black:background pixel (0), Blue: foreground pixel (1), Red:Don't care pixel (0 or 1)

Applied on a small skeleton-like test image (left), two branched-points  (junction) and four end-points are detected:

On a larger image, from the SIID images database:

Hit-or-miss transformation applied on the skeleton (up right), yields end-points (lower right). Since the same branched-points can be detected by different structuring elements, branched points appear of a variable color according to the detection occurence. Moreoever, branched-points can be connected yielding a branched domain of pixels.

ipython notebook

HMT-SE.ipynb

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.

 

Detecting end points or branched points from the skeleton of a single chromosome

The end-points or the branched-points from the skeleton of a single chromosome can be extracted using mathematical morphology operators implemented in mahotas.

Two images of metaphasic human chromosomes hybridized with a Cy3-telomeric PNA probe and  counterstained with DAPI, were used.

 

False color DAPI image

False color telomeres image

The location of a chromosome was identified, some preprocessing step were performed to segment the chromosome (hi-pass filtering followed by a simple thresholding).

The end-points, the branching points are extracted from the skeleton with a hit-or-miss operator implemented in mahotas, then overlayed on the original chromosome image.

The contour of the segmented telomeres are also overlayed:

Depending on the amount of high-pass filtering the chromosome skeleton fit the chromosome axis or the chromatids:

Skeletonization after "low" high-pass filtering

Skeletonization after a stronger high-pass filtering

In both cases, the skeleton is displayed in red, the branched-points in green, the end-points in blue, the chromosome contour in magenta, and finally the telomeres in yellow.
To run the code (python 2.7), the path to the images must be modified and mahotas and pymorph must be installed.
Download code

Getting the neighborhood of labelled particles in a graph

In the previous post, a labelled image, typically obtained by labelling connected components of a binary image, was used to compute a dictionary representing the neighborhood relationships between the particles.

A dictionary may be not the best data structure to represent such relationship. A graph structure is better suited and networkx provides a python implementation for that.

Starting from the following label image:

Background Grey level=0, first particle: grey level=1, ...last particle:grey level=5

We get two possible graphs. The first one takes the background (vertex 0) into account so the background is a neighbour of all the particles (there an edge between the vertex 0 and all the other vertex)

In this second graph, the background was removed:

A function  converts the dictionary to a networkx graph:

# -*- coding: utf-8 -*-
"""
Created on Mon Jun 25 10:26:01 2012

@author: Jean-Patrick Pommier
"""

import numpy as np
import networkx as nx
import mahotas as mh
import pylab as plb 

def makelabelarray():
    label = np.array([[0,1,4,4],
                     [1,0,0,2],
                     [3,0,2,2],
                     [0,2,0,5]])
    return label

def convertToGraph(dic, noBack=True):
    G = nx.Graph()
    G.add_nodes_from(dic.keys())
    for particle in dic.keys():
        list_touching_particles = dic[particle]
        # remove background
        if noBack:
            list_touching_particles.discard(0) 
        print 'v(',particle,')=',list_touching_particles
        for tp in list_touching_particles:
            G.add_edge(particle,tp)
    return G
    
def findneighborhoods(label,neighborhood):
    ''' given a labelled image, should return the adjacency list
        of particles, for a given neighborhood:
        
        neighborhood=np.array([0,1,0],[1,1,1],[0,1,0])
        
        The background (0), is kept as a particle neighbor 
        No fancy indexing
    '''
    #make the labels list
    labmax = label.max()
    #print labmax
    neighb_dic = {} # a dictionnary containing particle label as key and neighborhood
    for i in range(1,labmax+1):
        mask = (label ==i)
        #print mask
        dilated = mh.dilate(mask,neighborhood)
        neighbor = np.logical_and(dilated, np.logical_not(mask))
        #print neighbor
   #=======================================================================
        flatlab = np.ndarray.flatten(label)
        flatneighborhood = np.ndarray.flatten(neighbor)        
        flatneighbors = flatlab[flatneighborhood]
        flatneighbors.sort()
        #set is a trick so that each value of the neighborhoods is present only once
        neighb_dic[i] = set(flatneighbors)
        #print np.nonzero(flatneighbors)
    return neighb_dic
        
if __name__ == "__main__":
    a = makelabelarray()
    n = np.array([[1,1,1],[1,1,1],[1,1,1]])
    g = findneighborhoods(a, n)
    G = convertToGraph(g, noBack=True)

    plb.imshow(a,interpolation = 'nearest')
    plb.colorbar(ticks=[0,1,2,3,4])
    plb.show()

    nx.draw(G)
    plb.show()