# Nikroblog

Posts mostly around artifical intelligence and science.

Project maintained by Nikronic

# Ellipse Specific Fitting, RANSAC and Extracting Ellipse

This post consists of:

1. This step consists of following tasks:
1. Read This paper and summarize it
2. Convert MATLAB code in fig. 7 to python and extract ellipse parameters of `circle.bmp` and `ellipse.bmp` images.
3. Plot estimated ellipses using `cv2.ellipse()` method.
2. We want to find the parameters of ellipse using RANSAC algorithm. If only %40 of edges in the image belong to ellipse’s edges and we want to obtain the correct parameters with probability of 0.999, how many iterations are required?
3. Do these steps in this task:
1. Estimates parameters of ellipse using code in task 1 on `ellipse_noise.bmp` image
2. As there are points that do not belong to ellipse, RANSAC is better solution here. Implement RANSAC
3. Draw the output on `ellipse_noise.bmp` image
4. Set the probability of achieving correct parameters of ellipse to 0.99 and run algorithm for 10000 times. In how many of iterations, the estimated parameters are correct?

## 1 This step consists of following tasks:

1. Read This paper and summarize it
2. Convert MATLAB code in fig. 7 to python and extract ellipse parameters of circle.bmp and ellipse.bmp images.
3. Plot estimated ellipses using cv2.ellipse() method.

### 1.A Paper Summarization

The proposed method is ellipse specified which means no matter given data, a ellipse will be output. On top of that, it is computationally cheap and robust to noises. The major reason that this approach is robust and fast is that it uses least-square transformation.

First of all, they use a distance matrix with respect to ellipse equation which is called distance_matrix:

Ellipse Equation: Distance Matrix: Now, the parameter `a` is constrained using matrix called `C` which is 6x6 and all these constraints are linear or `C.dot(a) = 1`. But in this paper, constrained `a` is in the way that forces the fitted model to be ellipse. `4*a*c-b**2 = 1` is the equality constraint where `a.T.dot(C).dot(a) = 1`.

So `C` is: Based on what they have covered so far, the solution of the quadratically constrained minimization will be: Furthermore, this system can be written as below image where `lambda` is Lagrange multiplier and `S` is `D.T.dot(D)`: This system can be solved using generalized eigenvectors of `S.dot(a) = lambda*C.dot(a)`. In the end, if `(lambda, u)` solves `S.dot(a) = lambda*C.dot(a)`, we have: Then `a` can be obtained by `a = mu*u`.

### 1.B Direct Least Square of Fitting Ellipse Implementation and Ellipse of `circle.bmp` and `ellipse.bmp`

``````import numpy as np
import matplotlib.pyplot as plt
import cv2
%matplotlib inline
``````
``````def direct_least_square(x, y):
D = np.mat(np.vstack([x**2, x*y, y**2, x, y, np.ones(len(x))])).T
S = np.dot(D.T, D)
C = np.zeros((6, 6))
C[0, 2] = 2
C[1, 1] = -1
C[2, 0] = 2
Z = np.dot(np.linalg.inv(S), C)
eigen_value, eigen_vec = np.linalg.eig(Z)
eigen_value = eigen_value.reshape(1, -1)
pos_r, pos_c = np.where(eigen_value>0 & ~np.isinf(eigen_value))
a = eigen_vec[:, pos_c]
return a

def ellipse_center(a):
a = a.reshape(-1, 1)
b,c,d,f,g,a = a/2, a, a/2, a/2, a, a
num = b*b-a*c
x0=(c*d-b*f)/num
y0=(a*f-b*d)/num
return (int(y0[0, 0])+1, int(x0[0, 0])+1)

def ellipse_angle_of_rotation(a):
a = a.reshape(-1, 1)
b,c,d,f,g,a = a/2, a, a/2, a/2, a, a

def ellipse_axis_length(a):
a = a.reshape(-1, 1)
b,c,d,f,g,a = a/2, a, a/2, a/2, a, a
up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
res1=np.sqrt(up/down1)
res2=np.sqrt(up/down2)
return (int(res1[0,0]), int(res2[0, 0]))
``````
``````# read images

x_circle, y_circle = circle.nonzero()
x_ellipse, y_ellipse = ellipse.nonzero()

a_circle = direct_least_square(x_circle, y_circle)
a_ellipse = direct_least_square(x_ellipse, y_ellipse)
``````

### 1.C Plot Estimates

``````center = ellipse_center(a_circle)
axis = ellipse_axis_length(a_circle)
angle = ellipse_angle_of_rotation(a_circle)
start_angle = 0
end_angle = 360
color = 150
thickness = 1
plt.imshow(cv2.ellipse(circle, center, axis, angle, start_angle, end_angle, color, thickness))
``````
``````<matplotlib.image.AxesImage at 0x1749279ac50>
`````` ``````center = ellipse_center(a_ellipse)
axis = ellipse_axis_length(a_ellipse)
angle = ellipse_angle_of_rotation(a_ellipse)
start_angle = 0
end_angle = 360
color = 150
thickness = 1
plt.imshow(cv2.ellipse(ellipse, center, axis, angle, start_angle, end_angle, color, thickness))
``````
``````<matplotlib.image.AxesImage at 0x174927ee4a8>
`````` ## 2 How many Iterations for %40 Inlier Data With 0.999 Correct Estimation Probability?

``````w = 0.4
p = 0.999
# we need at least 6 points to estimates 6 parameters of the ellipse

k = np.log(1-p) / np.log(1-np.power(w, 6))
print('Number of needed iterations: {}'.format(int(np.ceil(k))))
``````
``````Number of needed iterations: 1684
``````

## 3 Do these steps in this task:

1. Estimates parameters of ellipse using code in task 1 on ellipse_noise.bmp image
2. As there are points that do not belong to ellipse, RANSAC is better solution here. Implement RANSAC
3. Draw the output on ellipse_noise.bmp image
4. Set the probability of achieving correct parameters of ellipse to 0.99 and run algorithm for 10000 times. In how many of iterations, the estimated parameters are correct?

### 3.A Estimate Ellipse on `ellipse_noise.bmp` Via Step 1 Code

``````# read images
x_ellipse_noise, y_ellipse_noise = ellipse_noise.nonzero()
a_ellipse_noise = direct_least_square(x_ellipse_noise, y_ellipse_noise)
``````
``````center = ellipse_center(a_ellipse_noise)
axis = ellipse_axis_length(a_ellipse_noise)
angle = ellipse_angle_of_rotation(a_ellipse_noise)
start_angle = 0
end_angle = 360
color = 150
thickness = 1
plt.imshow(cv2.ellipse(ellipse_noise, center, axis, angle, start_angle, end_angle, color, thickness))
``````
``````<matplotlib.image.AxesImage at 0x17492cb6390>
`````` ### 3.B Implement RANSAC for Ellipse

``````import random

def ransac(image, max_iter, threshold=5):
ellipse_noise = image
data = ellipse_noise
ics = []
best_ic = 0
best_model = None
xn, yn = data.nonzero()
nzero = [(x1,y1) for x1, y1 in zip(xn, yn)]
for epoch in range(max_iter):
ic = 0
sample = random.sample(nzero, 6)
a = direct_least_square(np.array([s for s in sample]), np.array([s for s in sample]))
for x, y in sample:
eq = np.mat(np.vstack([x**2, x*y, y**2, x, y, 1])).T
if np.abs(np.dot(eq, a.reshape(-1,1))) <= threshold:
ic += 1
ics.append(ic)
if ic > best_ic:
best_ic = ic
best_model = a
return a, ics

a, _ = ransac(ellipse_noise, 500, 5)
``````

### 3.C Draw the Estimated Ellipse Via RANSAC

``````center = ellipse_center(a)
axis = ellipse_axis_length(a)
angle = ellipse_angle_of_rotation(a)
start_angle = 0
end_angle = 360
color = 150
thickness = 1
plt.imshow(cv2.ellipse(ellipse_noise, center, axis, angle, start_angle, end_angle, color, thickness))
``````
``````<matplotlib.image.AxesImage at 0x1749492c2e8>
`````` ### 3.D If P=0.99, With 10000 Iteration, How Many Correct Estimations?

``````ellipse_noise = cv2.imread('ellipse_noise.bmp', 0)
a, ics = ransac(ellipse_noise, 1000, 5)
``````
``````center = ellipse_center(a)
axis = ellipse_axis_length(a)
angle = ellipse_angle_of_rotation(a)
start_angle = 0
end_angle = 360
color = 150
thickness = 1
plt.imshow(cv2.ellipse(ellipse_noise, center, axis, angle, start_angle, end_angle, color, thickness))
``````
``````<matplotlib.image.AxesImage at 0x17494a7e320>
`````` I do not know why sometimes my `direct_least_square` function, generates `0` parameters and sometimes `12`, so the could is unstable for high number of loops. If you have any fix for this, I would be so appreciated if you could handle it using a PR! Thank you ;-)

Project maintained by Nikronic Hosted on GitHub Pages — Theme by mattgraham