Thursday, April 10, 2014

In search of necklaces for a maximum area quadrilateral

The overlap of two chromosomes can be resolved by the four-point method provided that from these four points, a maximum area quadrilateral can be defined
Maximal quadrilateral overlaid on the overlap of two chromosomes

Previously that kind of quadrilateral was found from the 24 permutations of  four points:
4!=24 quadrilaterals formed by 4 points chosen randomly.

In these 24 permutations, some of them are derived from the others by circular permutations; these circular permutations define a necklace. The area of ​​quadrilaterals obtained by circular permutation of a necklace, is invariant.
The maximal area quadrilateral can be found in the set of the 6 necklaces.

To find a necklace, let's start from an ordered list of four points :
A, B, C, D
and swap two elements as (A,B), we get a necklace:
B, A, C, D
From four randomly chosen point, the following python function yields the six possible necklaces :
For example:

When run in ipython cell on sagemathcloud, the unique_necklaces function applied on four number took:

%timeit unique_necklaces(1,2,3,4)
100000 loops, best of 3: 15.2 µs per loop

Applied on four characters 'a', 'b', 'c', 'd'  , the function yields if the original 4-tuple is included (remove # in line 10):
print unique_necklaces('a','b','c','d')

[['a', 'b', 'c', 'd'], 
['b', 'a', 'c', 'd'], 
['c', 'b', 'a', 'd'], 
['d', 'b', 'c', 'a'], 
['a', 'c', 'b', 'd'], 
['a', 'd', 'c', 'b'], 
['a', 'b', 'd', 'c']]
Then a maximal area quadrilateral can be found beyond the six necklaces. As the number of quadrilaterals explored is smaller, the new implementation is faster than the previous one based on brute force:

     points=[(random.randint(0,20),random.randint(0,20)) for i in range(0,4)]
     %timeit maxAreaQuad(*points)
     %timeit maxAreaQuad_brutforce(*points)

     1000 loops, best of 3: 218 µs per loop
     1000 loops, best of 3: 748 µs per loop

It's funny to note that the area of the maximal area quadrilaterals seems to follow an asymetric distribution:

Ipython Notebook : here

Thursday, March 20, 2014

Combining five images into one rgb MFISH image using scikit-image implementation of HSV color space.

MFISH is a molecular cytogenetic technics using combination of fluorescent DNA probes used to recognize chromosomes and to detect some type of rearrangements.
In a previous post, five images corresponding to different fluorochromes were combined using matrix multiplication. The matrix coefficients were chosen empiracally and validated a-posteriori by the aspect of the rgb image produced.
Scikit-image is a python library implementing color space conversions for example between RGB and HSV. To display simultaneously the images of the five spectral components, or the five fluorochromes:

Five greyscaled images of the spectral components shown in false colors (cy55, cy5, SpGreen, SpOrange, TexasRed)
 The five greyscales images were first converted into 8bits rgb colors and their hue were modified with colors chosen prior combining them:
Greyscaled images converted into rgb images
Then the five images were added to produce one rgb MFISH image of type float64:

Let' produce the HSV histograms of the image:
Since the produced image lacks some contrast, the image exposition can be corrected such that the image background becomes black and the colors become brighter:

Now the histograms of the HSV chanels look like:
The linear stretch was also chosen based on the resulting rgb image.


The whole code is available in a ipython notebook.

Code on SageMathCloud:

The ipython notebook can be run on SageMathCloud. For some reasons, it's not possible to download images from github to sagemathcloud, so local images are used.

Images database:

Wednesday, February 5, 2014

Making image processing on SageMathCloud. Extraction of chromosomes from multispectral images.

SageMathCloud is a plateform to do mathematics, it also offers the possibility to use a ipython notebook online. Several python libraries are available: numpy, scipy, scikit-image (0.7, 0.93), mahotas 1.1, pandas ...
Once an account is opened (it's free), a project created, some multispectral images were uploaded in a directory of the project. In that project, an ipython notebook was created, this notebook can be run in a project directly in SageMathCloud.

Removing spurious pixels brought by ndimage.find_object():

The following model image  is a grey scaled image:

find.object() succeeds in extracting the letters:

An image of an extracted particle,for example the letter A, contain pixels belonging to other segmented (labelled) letters. The image of A contains pixels belonging to the letter r and to the letter B. Those spurious pixels have to be eliminated. find_object() provides the bounding box of the segmented particles in a list. The bouding box index seems to correspond to the particle label, this could be used to erase the spurious pixels:

A python function can be written:

It takes two images as arguments: a grey level image and a label image and returns a list of images containing the extracted particles from the grey level image. This implementation is much simpler and faster than a previous one .


Chromosomes extraction from a MFISH image:

After having uploaded images on cloud.sagemath, for example:

The five components (DAPI excepted) spectral images can be combined in one color image. There is an additional kar image resulting of the segmentation of the MFISH images :
5 spectral images combined into 3 to make a color (rgb) image (left). kar image (right)

Let's use the kar image to produce a label image so that each chromosome is defined by a single label (or grey level):
label image (right)

Now the extraction of individual chromosomes can be done, from the DAPI image:

or from stacked images:

The kar image identifies the chromosomes:

Let's make galleries of chromosomes

The extracted chromosome images have different size. Let's add borders to the images so that the images have all the same size :

Then let's build two galleries, one with MFISH images and one with inverse DAPI chromosomes, the code is visible in the notebook; it yields:
Left: MFISH chromosomes, Right:inverse DAPI

It remains to auto rotate the chromosome. A draft code, at the end of the ipython notebook  in the SageMathCloud project, tries to explore another method than the previously proposed. It remains also to arrange the chromosome according to their classification.

Download the ipython notebook

Tuesday, August 20, 2013

Moving chromosomes with kivy

A demo application available with kivy, was modified to move chromosomes:

The application depends on two files: and pictures.kv. Both files are in the same folder:
the particles folder contains the images to load.

The application should be modified so that the same zoom should be applied to all the chromosomes.

Update (22/08/2013) :

Independant zooming of the different chromosome images can be disabled.

Friday, August 16, 2013

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

Thursday, June 20, 2013

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


Wednesday, June 5, 2013

From binary image skeleton to networkx graph (work in progress ...)

To recognize objects in an image methods involving graphs have been developped as the shock graph method. These graphs are derived from image skeleton and then handled, compared. However, it seems taken for granted that a graph is available for work.
Here, a python code generating a networkx graph object from the skeleton of a binary image is proposed.

The image of the B letter was chosen as model. The image is generated with PIL:

Skeleton decomposition:

The skeleton was then computed by morphological thining, branched points were extracted with a hit-or-miss operator and labelled:

A bug somewhere (update fri 06, 2013):

For some shapes, as 'Y':

the algorithm fails at the end points linking step (see bellow). A way to cure the situation is to prune the initial skeleton of 'Y' once:
In the case of the letter 'A', the skeleton must be pruned three times for the algorithm to succeed:

Branched points neighborhood:

The end points can also be extracted and labelled from the initial skeleton, but end points from the edges will be needed.
The edges were obtained by difference between skeleton and branched points, there seven edges labelled from 1 to 7. The edges are pruned by one pixel (labels are 1,2,3,4,5,6,7).The edges end points were then extracted and labelled individually (1,2,3,...,14). The end points were also labelled according to their edge of origin ( 1,2,3,4,5,6,7):
For each branched point, its neighborhood was isolated to get the labels of the adjacent end points (right images):

Building the graph:

The graph was build by adding branched points first (vertex 1, 2, 3, 4), it assumes that at least one branched point exists, then  end points adjacent to the branched points are added (vertex 8,9,10 and 11,12,13 and 14,15,16): 
The remaining lonely end points are added to the graph (vertex 17, 18):
The end points are linked using the pruned edges, the edges are weighted according to their pixel size. The edges length here do not match to their weight:
Graph of the letter B, without initial skeleton pruning

Bug origin:

With the letter 'Y', with initial skeleton pruning, the algorithm yields the following graph:
Branched points are figured in red circles on the image of labelled edges and end points are figured by a white star.
In the case of the letter 'A', the algorithm yields :
The bug seems to come from too short edges as after three pruning steps removing the small barbs on the top of the letter, the algorithm succeeds in producing a graph

Download ipython notebook