how to plot ellipse given a general equation in R?

Manjunath picture Manjunath · Jan 24, 2017 · Viewed 10.9k times · Source

Ellipse general equation:

a * x ^ 2 + b * y ^ 2 + c * x * y + d * x + e * y + f = 0

enter image description here

Answer

Sandipan Dey picture Sandipan Dey · Jan 24, 2017

We can start from the parametric equation of an ellipse (the following one is from wikipedia), we need 5 parameters: the center (xc, yc) or (h,k) in another notation, axis lengths a, b and the angle between x axis and the major axis phi or tau in another notation.

enter image description here

xc <- 1 # center x_c or h
yc <- 2 # y_c or k
a <- 5 # major axis length
b <- 2 # minor axis length
phi <- pi/3 # angle of major axis with x axis phi or tau

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
plot(x,y,pch=19, col='blue')

enter image description here

Now if we want to start from the cartesian conic equation, it's a 2-step process.

  1. Convert the cartesian equation to the polar (parametric), form we can use the following equations to first obtain the 5 parameters using the 5 equations from the below figure (taken from http://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf, detailed math can be found there).

  2. Plot the ellipse, by using the parameters obtained, as shown above.

enter image description here

For step (1) we can use the following code (when we have known A,B,C,D,E,F):

M0 <- matrix(c(F,D/2,E/2, D/2, A, B/2, E/2, B/2, C), nrow=3, byrow=TRUE)
M <- matrix(c(A,B/2,B/2,C), nrow=2)
lambda <- eigen(M)$values
abs(lambda - A)
abs(lambda - C) 

# assuming abs(lambda[1] - A) < abs(lambda[1] - C), if not, swap lambda[1] and lambda[2] in the following equations:

a <- sqrt(-det(M0)/(det(M)*lambda[1]))  
b <- sqrt(-det(M0)/(det(M)*lambda[2]))
xc <- (B*E-2*C*D)/(4*A*C-B^2)
yc <- (B*D-2*A*E)/(4*A*C-B^2)
phi <- pi/2 - atan((A-C)/B)*2

For step (2) use the following code:

t <- seq(0, 2*pi, 0.01) 
x <- xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi)
y <- yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
plot(x,y,pch=19, col='blue')