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