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