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