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