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]

1 comment:

dip4fish said...

A more compact code from rroowwllaanndd:
http://stackoverflow.com/questions/6282462/converting-an-open-curve-to-a-list-of-ordered-pixels-a-python-test-code-with-num/6457666#6457666