Looking for bending points on the polygonal approximation of a contour

CurvedPointsOncontour.py

# -*- coding: utf-8 -*-
"""
Created on Tue Oct 23 08:06:37 2012

@author: Jean-Patrick
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import skimage.io as sio
import skimage
import skimage.viewer as viewer
from skimage.draw import ellipse
from skimage.measure import find_contours, approximate_polygon, \
    subdivide_polygon
import cmath
from math import pi

    
def addpixelborder(im):
    x,y = im.shape
    dy = np.zeros((y,))
    
    im = np.vstack((im,dy))
    im = np.vstack((dy,im))
    
    im = np.transpose(im)
    x,y = im.shape
    
    dy = np.zeros((y,))
    im = np.vstack((im,dy))
    im = np.vstack((dy,im))
    im = np.transpose(im)
    return im

def angle(p1,p2,p3):
    '''get the angle between the vectors p2p1,p2p3  with p a point
    '''
    z1=complex(p1[0],p1[1])
    z2=complex(p2[0],p2[1])
    z3=complex(p3[0],p3[1])
    v12=z1-z2
    v32=z3-z2
    print v12,v32
    angle=(180/pi)*cmath.phase(v32/v12)
    return angle
    
def angleBetweenSegments(contour):
    angles = []
    points = list(contour)
    points.pop()
    #print 'contour in def',len(contour)
    l = len(points)
    for i in range(len(points)+1):
        print 'i',i,' i-1',i-1,' (i+1)%l',(i+1)%l
        prv_point=points[(i-1)%l]
        mid_point=points[i % l]
        nex_point=points[(i+1)%l]
        
        angles.append(angle(prv_point,mid_point,nex_point))
    return angles




user=os.path.expanduser("~")
workdir=os.path.join(user,"QFISH","JPPAnimal","JPP52","11","DAPI","particles")
image='part25.png'
complete_path=os.path.join(workdir,image)
#print complete_path
im = sio.imread(complete_path)
im = addpixelborder(im)

fig, ax1= plt.subplots(ncols=1, figsize=(9, 4))


contours = find_contours(im, 0.01)
print type(contours)
print len(contours)
p = subdivide_polygon(contours[0], degree=5, preserve_ends=True)
#print p
coord = approximate_polygon(p, tolerance=1.5)
coord2 = approximate_polygon(p, tolerance=2)
print coord2[-1]
print coord2[0]
print coord2[1]

print 'coord',type(coord2[0])

ax1.imshow(im, interpolation='nearest', cmap=plt.cm.gray)

#plt.imshow(im, interpolation='nearest')
contour0=contours[0]
pf=contours[-1]
print 'contc',contour0[0]

angles = angleBetweenSegments(coord2)

table = np.array([coord2[0][0],coord2[0][1],angles[0]])
for i in range(len(angles[1:])):
    current =  np.array([coord2[i][0],coord2[i][1],angles[i]])
    table = np.vstack((table,current))

positive = table[table[:,2]>0]
negative = table[table[:,2]<0]

for n, contour in enumerate(contours):
    ax1.plot(contour[:, 1], contour[:, 0], linewidth=4)
    ax1.plot(coord[:, 1], coord[:, 0], linewidth=1, color='r')
    ax1.plot(coord2[:, 1], coord2[:, 0], linewidth=1, color='y')
ax1.scatter(contour0[0][1],contour0[0][0],s=150, c='m',marker='o')
ax1.scatter(pf[0][1],pf[0][0],s=150, c='y',marker='*')


for i in range(positive.shape[0]):
    ax1.scatter(positive[i][1], positive[i][0], s=200, c='r')
    
for i in range(negative.shape[0]):
    ax1.scatter(negative[i][1], negative[i][0], s=200, c='y')


plt.show()
#
#views = viewer.CollectionViewer((im))
#views.show()