how to plot ellipse given a general equation in R?

Question

Ellipse general equation:

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

Answer 1

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.

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

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.

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

Answer 2

The other answer shows you how to plot the ellipse, when you know both its centre and major axes. But they are not evident from the general ellipse equation. So here, I will start from scratch.

Omitting mathematical derivation, you need to solve for the centre from the following equation:

(oops: should be "generate v" not "generate u"; I can't fix it as the original LaTeX is now missing and I don't want to type again...)

Here is an R function to do this:

plot.ellipse <- function (a, b, c, d, e, f, n.points = 1000) {
  ## solve for centre
  A <- matrix(c(a, c / 2, c / 2, b), 2L)
  B <- c(-d / 2, -e / 2)
  mu <- solve(A, B)
  ## generate points on circle
  r <- sqrt(a * mu[1] ^ 2 + b * mu[2] ^ 2 + c * mu[1] * mu[2] - f)
  theta <- seq(0, 2 * pi, length = n.points)
  v <- rbind(r * cos(theta), r * sin(theta))
  ## transform for points on ellipse
  z <- backsolve(chol(A), v) + mu
  ## plot points
  plot(t(z), type = "l")
  }

Several remarks:

1. There are conditions for parameters a, b, ..., f in order to ensure that the equation is an ellipse rather than something else (say parabolic). So, do not pass in arbitrary parameter values to test. In fact, from the equation you can roughly see such requirement. For example, matrix A must be positive-definite, so a > 0 and det(A) > 0; also, r ^ 2 > 0.
2. I have used Cholesky factorization, as this is my favourite. However, the most beautiful result comes from Eigen decomposition. I will not pursue further on this part. If you are interested in it, read my another answer Obtain vertices of the ellipse on an ellipse covariance plot (created by car::ellipse). There are beautiful figures to illustrate the geometry of Cholesky factorization and Eigen decomposition.

Answer 3

You can use my package PlaneGeometry (on CRAN):

library(PlaneGeometry)

ell <- EllipseFromEquation(A = 4, B = 2, C = 3, D = -2, E = 7, F = 1)
box <- ell$boundingbox()
plot(NULL, asp = 1, xlim = box$x, ylim = box$y, xlab = NA, ylab = NA)
draw(ell, col = "yellow", border = "blue", lwd = 2)