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

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