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