## 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
```
Post a Comment