Fitted Ellipse and its major axis

[Adrián Galdrán]

Good night!

I have been trying to fit an ellipse to the border of a roughly elliptic 
region. The ellipse fits nicely, but whenever I try to plot the major axis, 
it fails. Of note, the approach below works well for phi \in [pi/4, 3pi/4], 
but not on other semi-cuadrants. I think the problem lies on the output 
of xc, yc, a, b, alpha = ellipse.params. Here, xc and yc are the 
coordinates of the center of the ellipse, a and b are the length of major 
and minor axis, and alpha is (I guess) the angle that the major axis makes 
against the horizontal axis. 

Does anyone have a hint for handling properly that angle, so that the code 
below works under any orientation of the ellipse? Many thanks!

from matplotlib import pyplot as plt

 

import numpy as np 

import skimage.io as io 

from skimage.measure import EllipseModel 

from skimage import feature 

from skimage.transform import rotate 

 

bin_im = io.imread(
'https://www.dropbox.com/s/fantqj8x4mbbsxf/bin_mask.png?dl=0') 

##################################################### 

#This routine fails for ellipses oriented in [0,pi/4] 

# bin_im = rotate(bin_im, 0, order=1) 

# bin_im = rotate(bin_im, 30, order=1) 


#It also seems to fail for [3*pi/4 ,pi] 

#bin_im = rotate(bin_im, 140, order=1) 
 

#However, it does not seem to fail for complementary 

#orientations belonging to [pi/4, 3*pi/4] - approximately 

# bin_im = rotate(bin_im, 70, order=1) 

#bin_im = rotate(bin_im, 100, order=1) 

##################################################### 



#Initial data processing 

fig, ax = plt.subplots(nrows=1, ncols=2) 

edges = feature.canny(bin_im) 

data_x, data_y = np.where( edges > 0.5 ) 

################################ 

#Subsample data points for speed 

indexes = np.arange(0,len(data_x),10) 

data_x = data_x[indexes] 

data_y = data_y[indexes] 

data = np.column_stack([data_x, data_y]) 

############################### 

#Plot original image and border 

ax[0].imshow(bin_im, cmap=plt.cm.gray) 

ax[1].imshow(edges, cmap=plt.cm.gray) 

plt.show() 

###################################### 

#Estimate ellipse parameters from data 

ellipse = EllipseModel() 

ellipse.estimate(data) 

xc, yc, a, b, alpha = ellipse.params 

######################################## 

#Build ellipse from estimated parameters 

#Rotation matrix 

R = np.array([ 

 [np.cos(alpha), -np.sin(alpha)], 

 [np.sin(alpha), np.cos(alpha)] 

]) 

a0, a1 = (0, 2*np.pi) 

angles = np.linspace(a0, a1, 300) 

X = np.vstack([ np.cos(angles) * a, np.sin(angles) * b]).T 

ell = np.dot(X, R.T) + (xc, yc) 

#################### 

#Plot fitted ellipse 

fig2, ax2 = plt.subplots(nrows=1, ncols=1) 

ax2.plot(bin_im.shape[1]-data[:, 1], data[:, 0], '.b', label='Input data') 

ax2.plot(bin_im.shape[1]-ell[:,1],ell[:,0],'r', label = 'Estimated data') 

ax2.legend(loc='upper left') 

# plt.show() 

#Equation of the lower half of the ellipse should be: 

#y = y_c + (x-x_c)*tan(alpha) 

major_axis_length = max(a, b) 

xx = np.linspace(xc-major_axis_length, xc+major_axis_length, 300) #as an 
aside question, how should i properly set this range?

yy = yc + (xx-xc)*np.tan(alpha) 

line = np.column_stack([xx, yy]) 

ax2.plot(bin_im.shape[1]-line[:, 1], line[:, 0], 'g') 

ax2.set_xlim([0, bin_im.shape[1]]) 

ax2.set_ylim([0,bin_im.shape[0]]) 

plt.show() 


 
#Build masked semi-ellipse 

h=bin_im.shape[0] 

w=bin_im.shape[1] 

xaxis = np.linspace(0, h-1, h) 

yaxis = np.linspace(0, w-1, w) 

yyy, xxx = np.meshgrid(yaxis, xaxis) 

mask = yyy > yc + (xxx-xc)*np.tan(alpha) 

new_bin_im = np.zeros_like(bin_im) 

new_bin_im[mask]=bin_im[mask] 

fig, ax3 = plt.subplots(nrows=1, ncols=1) 

ax3.imshow(new_bin_im,cmap=plt.cm.gray) 

plt.show()




-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scikit-image/attachments/20150512/7cb3252a/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bin_mask.png
Type: image/png
Size: 431 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/scikit-image/attachments/20150512/7cb3252a/attachment.png>