Monday, December 19, 2011

Telomeres quantification in a nucleus

Telomeres quantification can be done at low magnification in interphasic cells. A microscopic field was recorded with a CCD camera at low magnification (objective x25). Nuclei were counterstained with DAPI and the telomeres were detected with a Cy3-PNA probe.
From the field image, a nucleus is segmented and the DAPI and Cy3 spectral components are isolated. Even at low magnification, the telomeric spots can be segmented as follow:
Background is removed by top-hat filtering with a small structuring element (a disk of size=5). The top-hat image is then submitted to a hmin operator, leaving the peaks higher than 200 gray level units, then regional maxima are found.
The regional maxima are labelled to obtain markers for the watershed algorithm. Once thresholded, the top hat image is converted into distance map using euclidian or chess-board distance without significant difference.
Watershed doesn't succeed very well in separating spots belonging to the lower cluster of three spots. It may be interesting to use a high pass filtered image to build a  distance map.
The segmented spots can be overlaid to the original image (blue:telomeric spots, yellow:lines separating the spots after watershed). The quantification can be then performed on the original image (left) but the background must be estimated. The quantification can also be done on the Top-Hat filtered image where the background value is close to zero.
To quantify the spots, two parameters are measured: the spot area and the mean of the spot. The spot intensity is set to the product of the area by the mean gray level value (without taking the exposure time into account).
The area, mean and intensity values can be mapped onto an image or displayed as a histogram (from left to right: histogram of spots mean values, of spots area and histogram of spots intensity).

The following python  script relies on pymorph, mahotas and scipy. It could be tested with skimage.

# -*- coding: utf-8 -*-
"""
Created on Fri Dec 16 13:14:14 2011

@author: Jean-Patrick Pommier
"""

import scipy.ndimage as nd
import pylab as plb
import numpy as np

import mahotas as mh
import pymorph as pm

import skimage as sk
import skimage.morphology as mo

def E_distanceTransform(bIm):
    #from pythonvision.org
    dist = nd.distance_transform_edt(bIm)
    dist = dist.max() - dist
    dist -= dist.min()
    dist = dist/float(dist.ptp()) * 255
    dist = dist.astype(np.uint8)
    return dist
def Chess_distanceTransform(bIm):
    #from pythonvision.org
    dist = nd.distance_transform_cdt(bIm,metric='chessboard')
    dist = dist.max() - dist
    dist -= dist.min()
    dist = dist/float(dist.ptp()) * 255
    dist = dist.astype(np.uint8)
    return dist

def gray12_to8(im):
    i=0.062271062*im
    return pm.to_uint8(i)

path="/home/claire/Applications/ImagesTest/JPP50/16/"
CC="DAPI/"
PR="CY3/" 
im="1.TIF"

dapi=mh.imread(path+CC+im)
cy3=mh.imread(path+PR+im)
#==============================================================================
# Extracting known nuclei
#==============================================================================
nuc=dapi[215:264,1151:1198]
telo=cy3[215:264,1151:1198]
#==============================================================================
# Telomeres Segmentation with pymorph and ndimage
#==============================================================================
op=pm.open(telo,pm.sedisk(5))
top=telo-op
h=pm.hmin(top,200)
remax=pm.regmax(h)
blur=nd.gaussian_filter(h,3)
hi=h-1.0*blur
mask=np.copy(hi)
hi[mask<0]=0
remax=pm.edgeoff(remax)
markers,n=nd.label(remax)
print n," spots detected"
distE=E_distanceTransform(top>200)
distC=Chess_distanceTransform(top>200)
segE,L=pm.cwatershed(distE,markers,return_lines=True)
segC,L=pm.cwatershed(distC,markers,return_lines=True)

#==============================================================================
# Display segmentation on one nucleus
#==============================================================================
print 'display'
RedNuc=mh.bwperim((nuc>1000),n=4)
BlueSpot=mh.bwperim((segE>0),n=4)
GreenSpot=mh.bwperim((segC>0),n=4)


displayTelo=pm.overlay(gray12_to8(telo),red=RedNuc,blue=BlueSpot,yellow=L)
displayTeloC=pm.overlay(gray12_to8(top),red=RedNuc,green=GreenSpot)
#==============================================================================
# Quabtification
#==============================================================================
mean=pm.grain(top,segE,'mean','data')
print np.shape(mean)
area=pm.blob(segE,'area','data')
meanIm=pm.grain(top,segE,'mean','image')
areaIm=pm.blob(segE,'area','image')
print mean
plb.figure(1)
plb.subplot(141,frameon=False, xticks=[], yticks=[])
plb.title('raw 12bits')
plb.imshow(telo)
plb.subplot(142,frameon=False, xticks=[], yticks=[])
plb.title('top Hat')
plb.imshow(top)
plb.subplot(143,frameon=False, xticks=[], yticks=[])
plb.title('hmin')
plb.imshow(h)
plb.subplot(144,frameon=False, xticks=[], yticks=[])
plb.title('regional max')
plb.imshow(remax)

plb.figure(2)
plb.subplot(331,frameon=False, xticks=[], yticks=[])
plb.title('markers')
plb.imshow(markers)
plb.subplot(332,frameon=False, xticks=[], yticks=[])
plb.title('Thresholded Top-Hat')
plb.imshow(top>200)
plb.subplot(334,frameon=False, xticks=[], yticks=[])
plb.title('Euclid dist map ')
plb.imshow(distE)
plb.subplot(335,frameon=False, xticks=[], yticks=[])
plb.title('watershed seg')
plb.imshow(segE)
plb.subplot(333,frameon=False, xticks=[], yticks=[])
plb.title('TopH+hi pass')
plb.imshow(hi)
plb.subplot(337,frameon=False, xticks=[], yticks=[])
plb.title('Chess d map')
plb.imshow(distC)
plb.subplot(338,frameon=False, xticks=[], yticks=[])
plb.title('watershed Chess')
plb.imshow(segC)

plb.figure(3)
plb.subplot(121,frameon=False, xticks=[], yticks=[])
plb.title('telo+seg')
plb.imshow(displayTelo)
plb.subplot(122,frameon=False, xticks=[], yticks=[])
plb.title('telo+seg')
plb.imshow(displayTeloC)

plb.figure(4)
plb.subplot(231,frameon=False, xticks=[], yticks=[])
plb.title('Quant image=mean')
plb.imshow(meanIm)
plb.subplot(232,frameon=False, xticks=[], yticks=[])
plb.title('Quant img=area')
plb.imshow(areaIm)
plb.subplot(233,frameon=False, xticks=[], yticks=[])
plb.title('Quant img=mean*area')
plb.imshow(areaIm*meanIm)
plb.subplot(234)
plb.hist(mean)
plb.subplot(235)
plb.hist(area)
plb.subplot(236)
Q=areaIm*meanIm
plb.hist(Q.flatten()[Q.flatten()>0])
plb.show()

Friday, December 2, 2011

How to rotate faster chromosomes vertically (a little bit)

It is possible to orientate chromosomes vertically with a morphological  "rose des directions" using successive erosions. The larger the particle is, the longer the operation is. However, the iterations can be reduced by submitting not the particle itself, but a Top-Hat version of the image, or the particle skeleton or the particle contour yielding similar results with a reduced cost.

The first rose-des-directions was obtained with the original particle, I didn't compute the total amount of erosion, but finding the principal orientation around 120° required around 50 successive erosions, whereas only 16 erosions are necessary if the particle skeleton is analysed, 30 erosions with the Top-Hat filtered image or 20 with the contour.
With a more complex image, the "rose des directions" detects the orientation of the larger chromosome.

Thursday, November 3, 2011

Build Mamba 1.0 on Ubuntu 10.10

Mamba is a python mathematical morphology library. When dependencies are satisfied, installation is fast and simple. Once installed, from a python console:

>>> import mamba
>>> dir(mamba) 
['DEFAULT_GRID', 'EMPTY', 'FILLED', 'HEXAGONAL', 'MambaError', 'SQUARE', 'VERSION', '__builtins__', '__doc__', '__file__', '__name__', '__package__', '_always_show', '_blue_val', '_edge', '_green_val', '_grid', '_i', '_image_index', '_j', '_k', '_red_val', 'add', 'addConst', 'basinSegment', 'buildNeighbor', 'checkEmptiness', 'compare', 'computeDistance', 'computeMaxRange', 'computeRange', 'computeVolume', 'convert', 'convertByMask', 'copy', 'copyBitPlane', 'copyBytePlane', 'copyLine', 'cropCopy', 'diff', 'diffNeighbor', 'divConst', 'dualbuildNeighbor', 'extractFrame', 'generateSupMask', 'getDirections', 'getDisplayer', 'getHistogram', 'getImageCounter', 'getShowImages', 'gridNeighbors', 'hierarBuild', 'hierarDualBuild', 'hitOrMiss', 'imageMb', 'infFarNeighbor', 'infNeighbor', 'inverted_rainbow', 'label', 'logic', 'lookup', 'mambaCore', 'mbUtls', 'mul', 'mulConst', 'negate', 'os', 'patchwork', 'rainbow', 'raiseExceptionOnError', 'raiseWarning', 'rotateDirection', 'setDefaultGrid', 'setImageIndex', 'setMaxDisplay', 'setMinDisplay', 'setShowImages', 'shift', 'sub', 'subConst', 'supFarNeighbor', 'supNeighbor', 'threshold', 'tidyDisplays', 'transposeDirection', 'watershedSegment']

Thursday, October 20, 2011

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'

Monday, October 17, 2011

chromosome classification:workflow

The work-flow from raw DAPI image to two features scatter plot can be summarised as follow:
Yellow:to be written
Combining five metaphases, with two features (the area of the segmented particles and the ratio between the particle area and the area of convex hull of the particle) yields the following scatter plot where three clusters can be distinguished:
  • blue cluster: nuclei
  • red cluster: overlapping chromosomes
  • green cluster:single chromosomes
  • pink cluster:particles resulting from segmentation artefacts, some corresponds to nuclei touching the image border and should have been classified as nuclei.
As seen in the previous post the green and red clusters partially overlap. The graphic was build as follow with pandas (with some help):

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)
##
metalist=[0,1,2,3,4]
#metaphase=2
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',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')
pylab.show()

Thursday, October 13, 2011

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()

Wednesday, October 12, 2011

making a scatter plot with tags: playing with the iris data set

I need to display a scatterplot and to distinguish the points according to a category. I start with the iris dataset available in scikit-learn:
# -*- coding: utf-8 -*-
"""
"""
import pylab
from scikits.learn import datasets
iris = datasets.load_iris()
X=iris.data
Y=iris.target
print X
fig=pylab.figure()
ax = fig.add_subplot(111, aspect='equal')
print type(iris)
print iris.data.shape
print X[0:4,1]
print X[0:4,2]
print Y[0:4]
ax.scatter(X[:,0],X[:,1],c=Y[:])
pylab.show()

The Y column contains the labels stored as numerical values used to select the a color.

Friday, September 30, 2011

Chromosomes classification with a configurable pygame script :miniKar.py

 
Up to now miniKar doesn't recognize chromosomes. It is used to build a training set for :
  • supervised learning.
  • features selection
miniKar can be configured with a configuration file which specify the path to the chromosomes images and the categories of of the training set (but also the size and positions of the red areas, representing the categories).
#read by miniChromClassif.py
[Fluorochromes]
Counter_stain=DAPI
Probes=Cy3,Cy5,FITC
[Images Path]
work_dir=/home/claire/Applications/ImagesTest/
user=jp
slide=Jpp48
field_list=13,14
[ClassifierGUI]
screen_size=1024,768
categories_number=4
categories_list=single chr,touching chr,nuclei,dusts
box_size=240,380
#((col,li),(col+200+30,li),...)!!
box_position=(20,300),(268,300),(516,300),(764,300)
From left to right, the red areas represents
  • the single chromosomes.
  • the overlapping chromosomes
  • the nuclei
  • dusts (or segmentation artifacts)
By varying, the number of categories, it is already possible to perform by hand, a real karyotyping, even there is no tool to resolve overlapping chromosomes.
To use miniKar, three files must be in the same directory:
Suppose we want to classify the following particles into one of the four categories:
  1. single chromosome
  2. touching chromosomes
  3. nuclei
  4. dust
First download a chromosome dataset  then store it somewhere. On a ubuntu/linux the path to the images may be :
/home/user/cytogenetic_images/project/exampleslide/metaphase/DAPI/particles
Modify the configuration.cfg file to fit the path, for example:

work_dir=/home/user/cytogenetic_images
user=project
slide=exampleslide
field_list=1

Run minikar.py from a console by : python minikar01.py or from your favorite IDE (spyder, ninja-ide,...). Something like this should be visible:
The green background is uggly, but the small things are visible. After having classified the different chromosomes, one can have :
On pressing  the keyboard Esc,  minikar save a first file which associes the file name of a particle and its category and produce a second one in the csv format.
To classify the chromosomes from several  metaphases, the configuration.cfg file must be modified for each metaphase:
field_list=1
To classify the metaphase 13 in the folder exampleslide , modify the file as:
 field_list=13
If features are computed for each particle, then they can be associated to a category for further processing with a classifier. Here, chromosomes from 14 metaphases were classified, five features related to the convex hull were used.



Thursday, September 22, 2011

Let's try rotation

I am not too satisfied of the result. The interaction is not too convenient (key board+mouse), but for a first trial, that does the job. To rotate the chromosome, hit the r key of the keyboard to toggle between the drag-and-drop and the rotation modes.

Here is the script:

# -*- coding: utf-8 -*-
"""
Created on Mon Sep  5 15:44:42 2011

@author: jean-pat
"""
#!/usr/bin/env python
#
import os, pygame
import ConfigParser as CfgP
from pygame.locals import*

#def makeCfgFile():
#    config = CfgP.RawConfigParser()
#    config.add_section('shape classifier')
#    config.set('shape classifier', 'ref', '')
#    config.set('shape classifier', 'shape type', '')
#    return config
    
#def saveCfgFile(cfgfile,IKromList):
#    #building a config file
#    for s in IKromList:
#        print "saving config"
#        cref=s.ref
#        ctype=s.type
#        config.add_section(cref)
#        config.set(cref,'shape type', ctype)
#    config.write(open('shape.cfg','w'))
#    print "config saved"                
def readconfiguration():
    """
    read path to particles to classify:
        [Fluorochrome]
            Counter stain='DAPI'
            Probes=('Cy3','Cy5,'FITC')
        [Images Path]    
            work dir=...
            user=...
            slide=...
            field list=...
        [Classifier]
            screen size=(1024,768)
            category number=4
            category list=(single chromosome,touching chromosomes,nuclei,dust)
            box size=(200,100)
            first box=(5,650)
            box hspacing=30
 """
def load_image(name, colorkey=None):
    fullname=os.path.join("data", name)
    try:
        image=pygame.image.load(fullname)
    except pygame.error, message:
        print "Impossible de charger l'image:",name
        raise SystemExit, message
    image = image.convert()
    if colorkey is not None:
        if colorkey is -1:
            colorkey = image.get_at((0, 0))
        image.set_colorkey(colorkey, RLEACCEL)
    return image, image.get_rect()
    
                
class Ichrom(pygame.sprite.Sprite):
    def __init__(self,image,ref,initpos):
        pygame.sprite.Sprite.__init__(self)
        self.ref=ref#particule name
        self.type=""
        self.pos = initpos
        self.image,self.rect=image
        self.original=image
        self.rotation=0
        self.mode="translation"
        #print "self.rect",self.rect,"-self.image:",self.image
        self.button=(0,0,0)#mouse buttons not pressed
        #self.selected = 0
        self.mouseover=False
        self.focused=False
        self.rect.topleft=initpos
        print "init chrom at ",self.pos
    def rollover(self,rotationmode):
        """Test if the mouse fly over the chromosome 
        self.mouseover==True if mouse flying over the chrom,
        False if not"""
        mpos=pygame.mouse.get_pos()#mouseposition
        #test if mouse roll over the sprite
        if rotationmode==0:
            self.mode="translation"
        if rotationmode==1:
            self.mode="rotation"
            print pygame.mouse.get_rel()
        if self.rect.collidepoint(mpos):
            self.mouseover=True
        else:
            self.mouseover=False
            
    def update(self,classifiergroup,background):
        self.button=pygame.mouse.get_pressed()
        mpos = pygame.mouse.get_pos()
        self.selected=self.rect.collidepoint(mpos)
        #the mouse flies over a chromosome
        if (self.mouseover):
            if self.mode=="translation":
                #print "mouse pos:",mpos
                collision=pygame.sprite.spritecollide(self,classifiergroup,False)
                #collision should contains only one element            
                if len(collision)>0:
                    #print collision[0].category
                    self.type=collision[0].category
                    print "particle "+self.ref+" is classified to "+self.type
                if self.button==(1,0,0):     
                    pos = pygame.mouse.get_pos()
                    self.rect.center = pos
            elif self.mode=="rotation":
                #test mouse movement
                mrel=pygame.mouse.get_rel()
                if self.rotation>=360:
                    self.rotation=0
                    self.image=self.original
                elif mrel[0]>0:
                    self.rotation=self.rotation+10
                    self.image=pygame.transform.rotate(self.original[0],self.rotation)
                    self.rect=self.image.get_rect(center=self.rect.center)
                elif mrel[0]<0:
                    self.rotation=self.rotation-10
                    self.image=pygame.transform.rotate(self.original[0],self.rotation)
                    self.rect=self.image.get_rect(center=self.rect.center)
                
class Classifier(pygame.sprite.Sprite):
    '''When a chrom is moved is moved into a category '''
    def __init__(self,initpos,category):
        pygame.sprite.Sprite.__init__(self)
        self.category=category
        self.image = pygame.Surface((100,100))
        self.image.set_colorkey((0,0,0))
        self.image = self.image.convert_alpha()
        self.rect= self.image.get_rect()
        #pygame.draw.rect(screen, color, (x,y,width,height), thickness)
        pygame.draw.rect(self.image, (255,0,0,255), (0,0,100,100),1)
        self.rect.topleft= initpos
        
    def update(self): 
        pass        
        
def main():
    pygame.init()
    screen = pygame.display.set_mode((320,300))
    pygame.display.set_caption("Karyotyper")
    pygame.mouse.set_visible(True)
    
    background = pygame.Surface(screen.get_size())
    background = background.convert()
    background.fill((0, 0, 0))

    screen.blit(background,(0, 0))
    pygame.display.flip()
    
    i1=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part15.png", -1)
    i2=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part12.png", -1)
    chr1 = Ichrom(i1,"part15",(0,0))
    #chr1 = Krom(i1,(0,0))
    chr2=Ichrom(i2,"part12",(30,30))
    #chr2=Krom(i2,(30,30))
    categ1=Classifier((5,150),"single")
    categ2=Classifier((110,150),"overlapping")
    categ3=Classifier((215,150),"other stuff")
    
    allsprites = pygame.sprite.RenderPlain((chr1,chr2))
    allcategories=pygame.sprite.RenderPlain((categ1,categ2,categ3))
    clock = pygame.time.Clock()
    
    config = CfgP.RawConfigParser()
    spritelist=(chr1,chr2)
    rotationmode=0
    while 1:
        clock.tick(60)
        for event in pygame.event.get():
            if event.type == QUIT:
                return
            elif event.type == KEYDOWN and event.key == K_ESCAPE:
                #spritelist=(chr1,chr2)
                #building a config file
                #saveCfgFile(config,spritelist)
                for s in spritelist:
                    print "saving config"
                    cref=s.ref
                    ctype=s.type
                    config.add_section(cref)
                    config.set(cref,'shape type', ctype) 
                config.write(open('shape.cfg','w'))
                print "config saved"                
                return
            elif event.type == KEYDOWN and event.key == K_r:
                rotationmode=rotationmode+1
                if rotationmode==1:
                    print "rotation ON"
                if rotationmode==2:
                    print "rotation OFF"
                    rotationmode=0
            if event.type ==pygame.MOUSEBUTTONDOWN:
                #need to be modified to handle a list of chromosomes
                chr1.rollover(rotationmode)
                chr2.rollover(rotationmode)
                
        allsprites.update(allcategories,background)
        allcategories.update()
        ##
        ##Search which chromosome is moved
        ##into which category and classify  
        ##that chromosome in that category
#        collision=pygame.sprite.groupcollide(allcategories,allsprites,False,False,None)
#        for classified in collision.keys():
#            print classified
        screen.blit(background,(0,0))
        allsprites.draw(screen)
        allcategories.draw(screen)
        pygame.display.flip()
    
if __name__ == '__main__': main()

Tuesday, September 20, 2011

A mini chromosomes classifier with pygame

The following script handles a list of two chromosomes in order to set them in one of the two categories: "single chromosome", "overlapping chromosomes". The chromosomes are displayed in the upper part of the screen, there are three categories (red squares) available:"single chrom", "overlapping chrom" or "other stuff, ex nuclei". Just drag and drop the chromosomes in the red squares to classify them.

The script records the classification (press esc) in a configuration file located in the same directory. Here, there is two particles (part12.png or part15.png) obtained from a previous image segmentation. The classifier stores, in the file shape.cfg,the result as:

[part12]
shape type = overlapping

[part15]
shape type = single

# -*- coding: utf-8 -*-
"""
Created on Mon Sep  5 15:44:42 2011

@author: jean-pat
"""
#!/usr/bin/env python
#
import os, pygame
import ConfigParser as CfgP
from pygame.locals import*

#def makeCfgFile():
#    config = CfgP.RawConfigParser()
#    config.add_section('shape classifier')
#    config.set('shape classifier', 'ref', '')
#    config.set('shape classifier', 'shape type', '')
#    return config
    
def saveCfgFile(cfgfile,IKromList):
    for s in IKromList:
        print "saving config"
        cref=s.ref
        ctype=s.ref
        config.set('shape classifier', 'ref', cref)
        config.set('shape classifier', 'shape type', ctype)
        config.write(open('shape.cfg','w'))

def load_image(name, colorkey=None):
    fullname=os.path.join("data", name)
    try:
        image=pygame.image.load(fullname)
    except pygame.error, message:
        print "Impossible de charger l'image:",name
        raise SystemExit, message
    image = image.convert()
    if colorkey is not None:
        if colorkey is -1:
            colorkey = image.get_at((0, 0))
        image.set_colorkey(colorkey, RLEACCEL)
    return image, image.get_rect()
    
                
class Ichrom(pygame.sprite.Sprite):
    def __init__(self,image,ref,initpos):
        pygame.sprite.Sprite.__init__(self)
        self.ref=ref#particule name
        self.type=""
        self.pos = initpos
        self.image,self.rect=image
        #print "self.rect",self.rect,"-self.image:",self.image
        self.button=(0,0,0)#mouse buttons not pressed
        #self.selected = 0
        self.mouseover=False
        self.focused=False
        self.rect.topleft=initpos
        print "init chrom at ",self.pos
    def rollover(self):
        """Test if the mouse fly over the chromosome 
        self.mouseover==True if mouse flying over the chrom,
        False if not"""
        mpos=pygame.mouse.get_pos()#mouseposition
        #test if mouse roll over the sprite
        if self.rect.collidepoint(mpos):
            self.mouseover=True
        else:
            self.mouseover=False
    def update(self,classifiergroup,background):
        self.button=pygame.mouse.get_pressed()
        mpos = pygame.mouse.get_pos()
        self.selected=self.rect.collidepoint(mpos)
        #the mouse flies over a chromosome
        if (self.mouseover):
            #print "mouse pos:",mpos
            collision=pygame.sprite.spritecollide(self,classifiergroup,False)
            #collision should contains only one element            
            if len(collision)>0:
                #print collision[0].category
                self.type=collision[0].category
                print "particle "+self.ref+" is classified to "+self.type
            if self.button==(1,0,0):     
                pos = pygame.mouse.get_pos()
                self.rect.center = pos

class Classifier(pygame.sprite.Sprite):
    '''When a chrom is moved is moved into a category '''
    def __init__(self,initpos,category):
        pygame.sprite.Sprite.__init__(self)
        self.category=category
        self.image = pygame.Surface((100,100))
        self.image.set_colorkey((0,0,0))
        self.image = self.image.convert_alpha()
        self.rect= self.image.get_rect()
        #pygame.draw.rect(screen, color, (x,y,width,height), thickness)
        pygame.draw.rect(self.image, (255,0,0,255), (0,0,100,100),1)
        self.rect.topleft= initpos
        
    def update(self): 
        pass        
        #self.pos = (10,10)
        #pygame.draw.rect(self.image, (255,0,0,255), (10,50,100,100),2)
    def trainClassifier(self):
        """the particles moved inside the square are set into one
        category"""
        pass
        
def main():
    pygame.init()
    screen = pygame.display.set_mode((320,300))
    pygame.display.set_caption("Karyotyper")
    pygame.mouse.set_visible(True)
    
    background = pygame.Surface(screen.get_size())
    background = background.convert()
    background.fill((0, 0, 0))

    screen.blit(background,(0, 0))
    pygame.display.flip()
    
    i1=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part15.png", -1)
    i2=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part12.png", -1)
    chr1 = Ichrom(i1,"part15",(0,0))
    #chr1 = Krom(i1,(0,0))
    chr2=Ichrom(i2,"part12",(30,30))
    #chr2=Krom(i2,(30,30))
    categ1=Classifier((5,150),"single")
    categ2=Classifier((110,150),"overlapping")
    categ3=Classifier((215,150),"other stuff")
    
    allsprites = pygame.sprite.RenderPlain((chr1,chr2))
    allcategories=pygame.sprite.RenderPlain((categ1,categ2,categ3))
    clock = pygame.time.Clock()
    
    config = CfgP.RawConfigParser()
    while 1:
        clock.tick(60)
        for event in pygame.event.get():
            if event.type == QUIT:
                return
            elif event.type == KEYDOWN and event.key == K_ESCAPE:
                spritelist=(chr1,chr2)
                
                
                #building a config file
                for s in spritelist:
                    print "saving config"
                    cref=s.ref
                    ctype=s.type
                    config.add_section(cref)
                    config.set(cref,'shape type', ctype)
                    
                config.write(open('shape.cfg','w'))
                print "config saved"                
                return
            if event.type ==pygame.MOUSEBUTTONDOWN:
                #need to be modified to handle a list of chromosomes
                chr1.rollover()
                chr2.rollover()
                
        allsprites.update(allcategories,background)
        allcategories.update()
        ##
        ##Search which chromosome is moved
        ##into which category and classify  
        ##that chromosome in that category
#        collision=pygame.sprite.groupcollide(allcategories,allsprites,False,False,None)
#        for classified in collision.keys():
#            print classified
        screen.blit(background,(0,0))
        allsprites.draw(screen)
        allcategories.draw(screen)
        pygame.display.flip()
    
if __name__ == '__main__': main()

Monday, September 12, 2011

Drag and drop two chromosomes with pygame

Starting to understand how to move sprites with pygame, a small result:
The following code is still not a karyotyper (no classifier). It is a model to understand how to move two sprites.

# -*- coding: utf-8 -*-
"""
@author: jean-pat
"""
#!/usr/bin/env python
#
import os, pygame
from pygame.locals import*

def load_image(name, colorkey=None):
    fullname=os.path.join("data", name)
    try:
        image=pygame.image.load(fullname)
    except pygame.error, message:
        print "Impossible de charger l'image:",name
        raise SystemExit, message
    image = image.convert()
    if colorkey is not None:
        if colorkey is -1:
            colorkey = image.get_at((0, 0))
        image.set_colorkey(colorkey, RLEACCEL)
    return image, image.get_rect()

class Ichrom(pygame.sprite.Sprite):
    def __init__(self,image,initpos):
        pygame.sprite.Sprite.__init__(self)
        self.pos = initpos
        self.image, self.rect=image
        self.button=(0,0,0)#mouse buttons not pressed
        self.selected = 0
        self.mouseover=False
        self.focused=False
        print "init chrom at ",self.pos

    def rollover(self):
        """Test if the mouse fly over the chromosome 
        self.mouseover==True if mouse flying over the chrom,
        False if not"""
        mpos=pygame.mouse.get_pos()#mouseposition
        #test if mouse roll over the sprite
        if self.rect.collidepoint(mpos):
            self.mouseover=True
        else:
            self.mouseover=False
    def update(self):
        self.button=pygame.mouse.get_pressed()
        mpos = pygame.mouse.get_pos()
        self.selected=self.rect.collidepoint(mpos)
        #the mouse flies over a chromosome
        if (self.mouseover):
            #print "mouse pos:",mpos
            if self.button==(1,0,0):     # )  
                pos = pygame.mouse.get_pos()
                self.rect.midtop = pos

def main():
    pygame.init()
    screen = pygame.display.set_mode((400,300))
    pygame.display.set_caption("Karyotyper")
    pygame.mouse.set_visible(True)
    
    background = pygame.Surface(screen.get_size())
    background = background.convert()
    background.fill((0, 0, 0))

    screen.blit(background,(0, 0))
    pygame.display.flip()
    
    i1=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part15.png", -1)
    i2=load_image("/home/claire/Applications/ImagesTest/jp/Jpp48/13/DAPI/particules/part9.png", -1)
    chr1 = Ichrom(i1,(0,0))
    chr2=Ichrom(i2,(30,30))
    allsprites = pygame.sprite.RenderPlain((chr1,chr2))
    clock = pygame.time.Clock()

    while 1:
        clock.tick(60)
        for event in pygame.event.get():
            if event.type == QUIT:
                return
            elif event.type == KEYDOWN and event.key == K_ESCAPE:
                return
            if event.type ==pygame.MOUSEBUTTONDOWN:
                #need to be modified to handle a list of chromosomes
                chr1.rollover()
                chr2.rollover()
        allsprites.update()
        screen.blit(background,(0,0))
        allsprites.draw(screen)
        pygame.display.flip()
    
if __name__ == '__main__': main()

Monday, September 5, 2011

This is not a karyotyper

video
This is a trial with pygame, to move a chromosome with the mouse from a modified version of the chimp script. (Screencast with recordmydesktop + arista)
The next step is to handle several chromosomes, one by one.

Friday, July 1, 2011

Karyotyping with pygame?

In a cytogenetic application, it is necessary to interact with chromosomes in order to classify them that is to do a karyotype.The pygame library may be suitable for that. To do a karyotype, supposing that all the chromosomes of a metaphase are segmented, the following actions are required:
  • moving the chromosomes with the mouse (translations)
  • setting the chromosome orientation (rotations)
  • modifying the way the chromosomes are displayed (DAPI, inverse DAPI, color).
  • Setting the chromosome classification according to its location in the image.
The following code, directly adapted from a pygame tutorial, displays a moving chromosome (no mouse interaction) on a background image which may be used to classify the chromosomes:

Screenshot of the running script
# -*- coding: utf-8 -*-
"""
Éditeur de Spyder
"""
import scipy
import pygame as pyg
#============================================
class GameObject:
    def __init__(self, image, height, speed):
        self.speed = speed
        self.image = image
        self.pos = image.get_rect().move(0, height)
    def move(self):
        self.pos = self.pos.move(0, self.speed)
        if self.pos.right > 600:
            self.pos.left = 0
#==========================================
screen=pyg.display.set_mode((640,480))
chrom=pyg.image.load('/home/claire/Applications/
ImagesTest/jp/Jpp48/13/DAPI/particules/part3.png').convert()
backFile='/home/claire/Applications/ProjetPython/
sprite/backgroundClassifier.png'
background = pyg.image.load(backFile).convert()
screen.blit(background, (0, 0))
karyo=[]
chr1=GameObject(chrom, 1, 1)
karyo.append(chr1)
#===============================================
while 1:
    for event in pyg.event.get():
        if event.type in (pyg.QUIT, pyg.KEYDOWN):
            pyg.sys.exit()
        for chr1 in karyo:
            screen.blit(background, chr1.pos, chr1.pos)
        for chr1 in karyo:
            chr1.move()
            screen.blit(chr1.image, chr1.pos)
        pyg.display.update()
        pyg.time.delay(100)

Monday, June 20, 2011

Closed curve to set of ordered pixels: a python test

To transform a closed digital curve, remove one pixel, remember it and submit the opened curve to the following python code.
For example, the minimal closed curve:

gives:
closed curve
[[0 0 0 0 0]
 [0 0 1 0 0]
 [0 1 0 1 0]
 [0 1 0 1 0]
 [0 0 1 0 0]
 [0 0 0 0 0]]
index list
[array([1, 2]), 
array([2, 1]), 
array([3, 1]), 
array([4, 2]), 
array([3, 3]), 
array([2, 3]), 
array([1, 2])]
The fourth pixel:
[4 2]
Using functions, the code is a little bit cleaner, it relies only on numpy and scipy.ndimage for mathematical morphology (hit or miss), since only endpoints are needed here.If branched points are required, I will continue to use mahotas implementation of hit or miss operator because it handles the "don't care" pixels written with a "2" in the structuring element. Thus some of my scripts may not be run on windows, because mahotas has some installation problems due to freeimage on that plateform. I don't know if the scipy implementation of the hit or miss operator handle the "don't care" pixels.
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 20 09:08:38 2011

@author: Jean-Patrick
"""
import time
import numpy as np
from scipy import ndimage as nd

def opencurveToPixelsList(im):
    def scipyEndPoints(imcurve):
        se1=np.array([[0,0,0],[0,1,0],[0,1,0]])
        se2=np.array([[0,0,0],[0,1,0],[0,0,1]])
        EndPoints=np.zeros(imcurve.shape,dtype=np.uint)
        #perform hit & miss 
        for i in (0,3):
            hit1=nd.morphology.binary_hit_or_miss(imcurve,se1)
            se1=np.rot90(se1)
            EndPoints=EndPoints+hit1
            hit2=nd.morphology.binary_hit_or_miss(imcurve,se2)
            se2=np.rot90(se1)
            EndPoints=EndPoints+hit2
        return EndPoints
    #total number of pixels in the curve
    pixel=np.array([-1,-1])
    pixN=np.sum(im==1)
    pixelsList=[]
    #get end point by H&Miss op by ndimage
    #print scipyEndPoints(im)
    #get end points by hit&miss operator provided by mahotas
    #ep=bep.get_endpoints(im)
    ep=scipyEndPoints(im)
    lab,n=nd.label(ep)
    #where the first endpoint is
    first_indices=np.where(lab==1)
    pixel[0]=first_indices[0][0]
    pixel[1]=first_indices[1][0]
    pixelsList.append(np.copy(pixel))
    #print "first point",first_indices," vector",walkingPixel
    firstEndPoint=np.uint8(lab==1)
    #keep an image of the second end point
    lastEndPoint=np.uint8(lab==2)
    last_Indices=np.where(lastEndPoint==1)
    #print "last point",last_Indices
    ######################################################
    ## walk on curve with a 3x3 neighborhood
    ######################################################
    #Init
    current_curve=np.copy(im)
    current_point=np.copy(firstEndPoint)
    ###start to walk on curve
    #remove the second last point
    #current_curve=current_curve-lastEndPoint##too heavy?, try indices
    current_curve[last_Indices[0][0],last_Indices[1][0]]=0
    c_point_ind=np.where(current_point==1)
    #print c_point_ind
    li=c_point_ind[0][0]
    col=c_point_ind[1][0]
    for i in range (0,pixN-2):    
        #3x3 neighborhood arround the endpoint
        neighbor=current_curve[li-1:li+2,col-1:col+2]
        neighbor[1,1]=0
        #################
        #Can only handle a curve
        #such np.where(neihgbor==1) must gives only one pixel
        #############
        nextPointIndices=np.where(neighbor==1)##vectO'M=nextPointIndices
        ##remove the first point from the curve
        current_curve[li,col]=0
        pixN=pixN-1
        #print current_curve
        ##compute nextPoint indices in the original base vectOM=OO'+O'M
        li=(li-1)+nextPointIndices[0][0]
        col=(col-1)+nextPointIndices[1][0]
        pixel[0]=li
        pixel[1]=col
        pixelsList.append(np.copy(pixel))
    #don't forget the last pixel
    pixel[0]=last_Indices[0][0]
    pixel[1]=last_Indices[1][0]
    pixelsList.append(np.copy(pixel))
    return pixelsList
############################
def closedcurveToPixelsList(im):
    '''
    given a closed curve (no test),
    not touching the image border(no test) as:
        [0,0,0,0,0],
        [0,0,1,0,0],
        [0,1,0,1,0],
        [0,1,0,1,0],
        [0,0,1,0,0],
        [0,0,0,0,0]
    The function returns an array list of ordered 
    points along the curve:
        [1,2],[2,3],[3,3],[4,2],[3,1],[2,1],[1,2]
    '''
    firstpoint=np.copy(np.where(im==1))
    print np.where(im==1)[1][0]
    print "first point",firstpoint
    openedcurve=np.copy(im)
    openedcurve[firstpoint[0][0],firstpoint[1][0]]=0
    print "now it's open"
    print openedcurve
    openlist=opencurveToPixelsList(openedcurve)
    #add the first point of the closed curved at the 
    #begiining and the end of the list:
    pixel=np.array([-1,-1])
    pixel[0]=firstpoint[0][0]
    pixel[1]=firstpoint[1][0]
    openlist.append(np.copy(pixel))
    openlist.insert(0,np.copy(pixel))
    return openlist
#tests
im=np.array([[0,0,0,0,0,0,0],
            [0,1,0,0,1,0,0],
            [0,0,1,0,0,1,0],
            [0,0,0,1,1,0,0],
            [0,0,0,0,0,0,0]])
liste=opencurveToPixelsList(im)
closedcurve=np.array([
        [0,0,0,0,0],
        [0,0,1,0,0],
        [0,1,0,1,0],
        [0,1,0,1,0],
        [0,0,1,0,0],
        [0,0,0,0,0]])
liste2=closedcurveToPixelsList(closedcurve)
print liste
print liste[0]
print "closed curve"
print closedcurve
print liste2
print liste2[3]

Friday, June 17, 2011

open curve to ordered pixels: a second test

The previous code is modified :
  • to restrict the search area in a 3x3 domain around the end points of the curve
  • to depends only on numpy and scipy
With a model curve looking like:
The script returns an ordered list of pixels:
listpix [
array([1, 1]),
array([2, 2]),
array([3, 3]),
array([3, 4]),
array([2, 5]),
array([1, 4])]
The version seems to be faster than the previous one:
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 16 23:43:15 2011

@author: Jean-Patrick Pommier
"""
import time
import numpy as np
from scipy import ndimage as nd
#import BranchedEndPoints as bep

def scipyEndPoints(imcurve):
    se1=np.array([[0,0,0],[0,1,0],[0,1,0]])
    se2=np.array([[0,0,0],[0,1,0],[0,0,1]])
    EndPoints=np.zeros(imcurve.shape,dtype=np.uint)
    #perform hit & miss 
    for i in (0,3):
        hit1=nd.morphology.binary_hit_or_miss(imcurve,se1)
        se1=np.rot90(se1)
        EndPoints=EndPoints+hit1
        hit2=nd.morphology.binary_hit_or_miss(imcurve,se2)
        se2=np.rot90(se1)
        EndPoints=EndPoints+hit2
    return EndPoints
##original image
im=np.array([[0,0,0,0,0,0,0],
            [0,1,0,0,1,0,0],
            [0,0,1,0,0,1,0],
            [0,0,0,1,1,0,0],
            [0,0,0,0,0,0,0]])
print im

#total number of pixels in the curve
pixel=np.array([-1,-1])
pixN=np.sum(im==1)
pixelsList=[]
#get end point by H&Miss op by ndimage
#print scipyEndPoints(im)
#get end points by hit&miss operator provided by mahotas
#ep=bep.get_endpoints(im)
ep=scipyEndPoints(im)
lab,n=nd.label(ep)
#where the first endpoint is
first_indices=np.where(lab==1)
pixel[0]=first_indices[0][0]
pixel[1]=first_indices[1][0]
pixelsList.append(np.copy(pixel))
#print "first point",first_indices," vector",walkingPixel
firstEndPoint=np.uint8(lab==1)
#keep an image of the second end point
lastEndPoint=np.uint8(lab==2)
last_Indices=np.where(lastEndPoint==1)
#print "last point",last_Indices
######################################################
## walk on curve with a 3x3 neighborhood
######################################################
#Init
current_curve=np.copy(im)
current_point=np.copy(firstEndPoint)
###start to walk on curve
#remove the second last point
#current_curve=current_curve-lastEndPoint##too heavy?, try indices
current_curve[last_Indices[0][0],last_Indices[1][0]]=0
c_point_ind=np.where(current_point==1)
#print c_point_ind
li=c_point_ind[0][0]
col=c_point_ind[1][0]
for i in range (0,pixN-2):    
    #3x3 neighborhood arround the endpoint
    neighbor=current_curve[li-1:li+2,col-1:col+2]
    neighbor[1,1]=0
    #################
    #Can only handle a curve
    #such np.where(neihgbor==1) must gives only one pixel
    #############
    nextPointIndices=np.where(neighbor==1)##vectO'M=nextPointIndices
    ##remove the first point from the curve
    current_curve[li,col]=0
    pixN=pixN-1
    #print current_curve
    ##compute nextPoint indices in the original base vectOM=OO'+O'M
    li=(li-1)+nextPointIndices[0][0]
    col=(col-1)+nextPointIndices[1][0]
    pixel[0]=li
    pixel[1]=col
    pixelsList.append(np.copy(pixel))
#don't forget the last pixel
pixel[0]=last_Indices[0][0]
pixel[1]=last_Indices[1][0]
pixelsList.append(np.copy(pixel))
print "listpix",pixelsList

Wednesday, June 8, 2011

Converting an open curve to a list of ordered pixels: a python test code

In some conditions, a particle skeleton is an open curve with no branched point, according to a 4-connected neighborhood:
Test curve
Getting the pixels indices according to their order on the curve may be usefull (computing direction, bending,fitting...).
The following SkelToPixel_2 script, depends on two other scripts:BranchedEndPoints and StructElem. The scripts must be stored in the same folder.
The idea is simple, the script takes the curve at one end, walk on it searching a neighbor pixel up to the ...antepenultimate pixel (that works).Opencv yields the coordinates of pixels belonging to the contour of a particle, i.e. a closed curve. I don't know how to, but I suppose it could be possible to do the same thing with opencv for an open curve.
# -*- coding: utf-8 -*-
"""
Created on Tue Jun  7 10:35:17 2011
                      SkelToPixel_2
A code to convert an open curve stored in a numpy array
into a list of neihbor pixels

@author: jean-Patrick Pommier
http://dip4fish.blogspot.com
"""

import numpy as np
from scipy import ndimage as nd
import BranchedEndPoints as bep

im=np.array([[0,0,0,0,0,0,0],
            [0,1,0,0,0,0,0],
            [0,0,1,1,0,0,0],
            [0,0,0,0,1,1,0],
            [0,0,0,0,0,0,0]]) 
#total number of pixels in the curve
pixN=np.sum(im==1)
pixelsList=[]
#get end points by hit&miss operator
ep=bep.get_endpoints(im)
lab,n=nd.label(ep)
#where the first endpoint is
first_indices=np.where(lab==1)
firstEndPoint=np.uint8(lab==1)
#keep an image of the second end point
lastEndPoint=np.uint8(lab==2)
pixelsList.append(first_indices)
current_curve=np.copy(im)
current_point=firstEndPoint
print current_curve
#start to walk on the curve
while pixN<>2:
    current_curve=current_curve-current_point
    pixN=pixN-1
    #this is not very smart, using a 3x3 neighborhood
    #arround the current point would reduce the research
    current_point=bep.get_endpoints(current_curve)-lastEndPoint
    pixelsList.append(np.where(current_point==1))
#add the last pixel
pixelsList.append(np.where(lastEndPoint==1))
for i in range(0,len(pixelsList)):
    print "i=",i,pixelsList[i]
Run from spyder, the script output is:
[[0 0 0 0 0 0 0]
 [0 1 0 0 0 0 0]
 [0 0 1 1 0 0 0]
 [0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0]]
i= 0 (array([1]), array([1]))
i= 1 (array([2]), array([2]))
i= 2 (array([2]), array([3]))
i= 3 (array([3]), array([4]))
i= 4 (array([3]), array([5]))
The script doesn't yield directly pixels indices.
Feel free to propose a nicer, faster code.

The result with another curve
And the output console:
[[0 0 0 0 0 0 0]
 [0 1 0 0 1 0 0]
 [0 0 1 0 0 1 0]
 [0 0 0 1 1 0 0]
 [0 0 0 0 0 0 0]]
i= 0 (array([1]), array([1]))
i= 1 (array([2]), array([2]))
i= 2 (array([3]), array([3]))
i= 3 (array([3]), array([4]))
i= 4 (array([2]), array([5]))
i= 5 (array([1]), array([4]))

Monday, June 6, 2011

A pig chromosomes gallery: a modified makeMosaic function

I fixed a bug in the makeMosaic function in the script used to build a chromosomes gallery. Instead of removing, too small or to large particles from a list, now the makeMosaic function build a new list of filtered particles (chromosomes). The script was modified to display chromosomes in inverse DAPI with
pylab.set_cmap(pylab.cm.Greys)
Inverse DAPI pig chromosomes
Original pig DAPI stained metaphase

Tuesday, May 31, 2011

A minimal use of opencv: toward chromosome shape analysis.

In an attempt to analyse the chromosomes shape, I need to access to the coordinates of points belonging to the contour or to the skeleton of chromosomes. Opencv is a rich vision library, usable from python. Here, I wrote a minimal script to access the pixels of the chromosome contour:
part15.png:image used to test opencv

import cv
definition of some colors
_red =  (0, 0, 255, 0);
_green =  (0, 255, 0, 0);
_white = cv.RealScalar (255)
_black = cv.RealScalar (0)
# create window and display the original picture in it
path="/home/claire/Applications/ProjetPython/test opencv/part15.png"
image = cv.LoadImage(path,cv.CV_LOAD_IMAGE_GRAYSCALE)
binary=image>0
cv.NamedWindow ("image", 1)
cv.ShowImage ("image", image)
# create the storage area
storage = cv.CreateMemStorage (0)
# find the contours
contours = cv.FindContours(image,
                               storage,
                               cv.CV_RETR_TREE,
                               cv.CV_CHAIN_APPROX_SIMPLE,
                               (0,0))
cv.DrawContours(image,contours,_red,_black,max_level=2)
print len(contours)
print contours[0],contours[1],contours[61]
cv.WaitKey (0)

  • The first thing I noticed is that without cv.WaitKey(0), the image is not visible, it is displayed and erased immediately.
  • The object contours is a list of pixels coordinates, I didn't check yet if it is an ordered list of pixels.
My feeling about opencv is that it is a world itself, may be not mixing very well with other library.

Monday, May 23, 2011

Making a MFISH image

The first trial with numpy was not very successful. After code simplification, I get a result:
RGB image combining five monospectral images
The code is:
# -*- coding: utf-8 -*-
"""
Created on Mon May 23 14:07:40 2011

@author: jean-pat
"""

import numpy as np
import os
import pylab
#from scipy import linalg as la
#from scipy import ndimage as nd
import scipy as sp
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=sp.misc.imread(complete_path)

R=0.8*i5+0.15*i4+0.05*i3+0.00*i2+0.00*i1
V=0.00*i5+0.05*i4+0.05*i3+0.80*i2+0.10*i1
B=0.00*i5+0.00*i4+0.05*i3+0.15*i2+0.80*i1
shape=i5.shape
rgb = np.zeros((shape[0],shape[1],3),dtype=np.uint8)
rgb[:,:,0]=np.copy(R)
rgb[:,:,1]=np.copy(V)
rgb[:,:,2]=np.copy(B)
pylab.subplot(111,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb)
pylab.show()

I received new advices today
# -*- coding: utf-8 -*-
"""
Created on Mon May 23 14:07:40 2011

@author: jean-pat
"""

import numpy as np
import os
import pylab
#from scipy import linalg as la
#from scipy import ndimage as nd
import scipy as sp
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=sp.misc.imread(complete_path)
######################
multi_spectr = np.dstack([i1,i2,i3,i4,i5])
conv_matrix = np.array([
[0.8,0.15,0.05,0,0],
[0,0.05,0.05,0.8,0.1],
[0,0,0.05,0.15,0.8]
])
 
shape = multi_spectr.shape
multi_spectr = multi_spectr.reshape(shape[0]*shape[1],shape[2])
rgb = np.dot(multi_spectr,conv_matrix.T).reshape(shape[:2]+(3,))
rgb8=np.uint8(rgb)
########################
#R=0.8*i5+0.15*i4+0.05*i3+0.00*i2+0.00*i1
#V=0.00*i5+0.05*i4+0.05*i3+0.80*i2+0.10*i1
#B=0.00*i5+0.00*i4+0.05*i3+0.15*i2+0.80*i1
#shape=i5.shape
#rgb = np.zeros((shape[0],shape[1],3),dtype=np.uint8)
#rgb[:,:,0]=np.copy(R)
#rgb[:,:,1]=np.copy(V)
#rgb[:,:,2]=np.copy(B)
pylab.subplot(111,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb8)
pylab.show()

That yields:

Thursday, May 19, 2011

MFISH image to RGB with image5D

Wade Schwartzkopf published a MFISH dataset  in 2000. Image5D (a plugin for imagej) can display five different fluorochromes simultaneously:
Left:raw png images. Top right: RGB image done with Image5D. Bottom right: "karyotype" image with false colors. 

Far red
Red
Spectrum gold
Spectrum green
Spectrum aqua
Dapi
After chromosome classification, a karyotype image can be produced. This is a grayscaled image associated to a  LUT, to be displayed with false colors:
Image5D can export a rgb image from the five channels image:
RGB image generated after background substraction of the five raw images

Tuesday, May 17, 2011

Highlighting Interstitial Telomeric Spots in interphasic cells (HITS-FISH)

Telomeric fusions underscore telomeric binding proteins disfunctions. Fusions occur in senescent cells between short telomeres, they are detected in metaphasic cells by dicentric chromosomes most generally with no telomeric spots at the fusion point.
In some SV40 transformed fibroblast cell line, interstitial telomeric spots (ITS) are as bright as terminal telomeric spots, this indicates some uncoupling between telomere lenght and their capacity to block telomeric fusions. 
ITS observed in a precrisis human SV40 transformed cell line (Ducray et al., oncogene 1999)

From metaphasic chromosomes observation, fusions retaining visible telomeric spots are rare events but may be more frequent in interphasic cells.
ITS from a SV40 transformed cell line

Bal31 in-situ:

Several nucleases or DNA polymerases  were used in assays requiring fixed material (chromosome banding+restriction endonuclease, in-situ nick translation, fluorescent labelling with TUNEL, in-situ PCR). The exonuclease Bal31 removes totally terminal telomeric sequences of naked DNA after some hours of digestion (R C Allshire, M Dempster, and N D Hastie, NAR 1989; de Lange et al, MCB 1990) . In-situ,  the terminal telomeric sequences, should be removed by Bal31 digestion.  If telomeric spots can be detected in nuclei after Bal31 digestion, those spots should be  ITS. Monitoring ITS may be a tool of interest  to probe telesome functions at the scale of unique cell, or to monitor the state of a tissue.

Pig cells as model cells?

On human metaphasic chromosomes, a bright ITS indicates an  abnormal cell but for cheaper DAPI stained chromosomes rearrangements too. Normal pig karyotype  shows interstitial telomeric spots:
Normal diploid pig metaphase with a chromosomes pair having ITS.
After QFISH on interphasic pig nuclei, there are two ITS of same magnitude than terminal spots, these two pig ITS remain hidden in the forest of terminal telomeric spots:
  
Pig nuclei (25x magnification)
Some conditions of Bal31 digestion should retain these two ITS and remove totally terminal telomeric spots yielding the possibility to observe telomeric fusions in nuclei.

Assay:
Tuning the conditions of Bal31 digestion is a prerequisite step. Supposing that the exonuclease Bal31 can remove terminal sequence from fixed material, when this enzyme is usually used to digest purified DNA in solution. 

If pig cells are not available, performing Bal31 digestion on control human cells such lymphocytes should tell if Bal31 can remove telomeres in-situ.
From conventionnal cytogenetic preparation of pig cells (fibroblasts,lymphocytes), digestion time, numbers of Bal31 units could be set up as follow : 
Possible protocol for HITS FISH

The protocol modified from Lansdorp et al. (Hum Mol Genet. 1996 May;5(5):685-91.) could be also adapted to liquid hybridization for flow FISH. Under optimal conditions pig cells could be used as internal control for Bal 31 digestion by mixing human and pig cells in the same cytogenetic preparation (9:1 ratio).

Human cells can be then distinguished from pig cells with an additionnal oligo probe. In the following example (25x magnification), simultaneous hybridization with Cy3-PNA-(CCCTAA)3 (displayed in green) and Cy5-DNA (displayed in red) oligo probe was performed (unpublished data):

Candidate cells to test the assay, may be SV40 precrisis cells with highly rearranged chromosomes and numerous dicentrics, or better cells with controlled mutations.