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

Solution 1

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.

Solution 2

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

Solution 3

You can use my package PlaneGeometry (soon on CRAN, hopefully):

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)

Comments

[curious_cat]

Great answer! Question: How come the original question had an 'f' parameter but your solution does not? Does it drop out?

[curious_cat]

My bad. I only looked at the first matrix equation.

[curious_cat]

So if I add a check on positive definiteness then that ensures that whatever arbitrary parameters a user enters those are ok for an ellipse?

[Zheyuan Li]

@curious_cat These must hold: a > 0; 4ab - c^2 > 0; also, r^2 > 0. This is a working example: plot.ellipse(1, 3, 2, 12, -5, 50)

[java_xof]

three remarks regarding the code itself: y <- yc + a * cos(t) * sin(phi) + b * sin(t) * cos(phi) and phi <- (pi/2 - atan((A-C)/B))/2 and xc <- (B*E-2*C*D)/(4*A*C-B^2)

[Adam]

error : the difference between equation ( image at the top of the post) and the code for y

[Adam]

removed error from code for y in 2 places

[Karl Forner]

for step2 I believe that you have a typo, it should be: y <- yc + a*cos(t)*sin(phi) + b*sin(t)*cos(phi)