Contour : looking for negatively curved domains

A contour was previously explored with the convexity defects  algorithm implemented in opencv. Some highly curved domains  failed to be detected because of the contour shape.

Curved domains missed by convexity defects detection

Using polygonal approximation imlemented in skimage 0.8 dev, positively curved domains as well as negatively curved domains on a contour can be sorted:

Contour of a chromosomes cluster (blue), polygonal approximation (yellow). Red points: positive curvature. Yellow points: negative curvature.

From a segmented particle, a polygonal approximation of the contour is extracted.

Points of the polygon are used to compute an angle between two consecutive vectors (Here using the complex affix of vectors):

Positive angles correspond to convex domains and negative angles to concave domains.

The code has to be modified to filter angles according to a given range. 

The code can be modified as follow to filter angles corresponding to curved domains:

• positive = table[(135>table[:,2]) & (table[:,2]>45)]

• negative = table[(-45>table[:,2]) & (table[:,2]>-135)]

Curved domains on polygonal approximation of a chromosomes cluster contour.

Download code

	# -*- 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()

view raw CurvedPointsOncontour.py hosted with ❤ by GitHub