R - calculating overlap between two ellipses

Question

First things first, yes, there have been similar questions asked here previously. However, they all stem from PCA or from some specialized packages; it's visible in the data that these people provided in the question that this data doesn't match the case that I have and I cannot use it.

I have a dataset literally with the centroids, major and minor axes, and angles of the ellipse. Here is a very minimal example:

data <- data.frame(x0 = c(0, 0), y0 = c(0, 0), a = c(5, 3), b = c(10, 20), angle = c(45, 35), Ellipse = c("Ell1", "Ell2"))

The dataframe:

  x0 y0 a  b angle Ellipse
1  0  0 5 10    45    Ell1
2  0  0 3 20    35    Ell2

I also display it visually, just for didactic purposes (I don't need to plot the overlapping area):

library(ggplot2)
library(ggforce)

ggplot(data, aes(x0 = x0, 
                 y0 = y0,
                 a = a, 
                 b = b, 
                 angle = angle, 
                 color = Ellipse)) + 
  geom_ellipse()

Given this setup, how can I calculate the overlapping area?

Answer 1

The math for an exact solution is very difficult here. I think the easiest way to do this in R is to create your ellipses as polygons and use the sf library to get their intersection. It is then very easy to get its area.

The most difficult part is creating your ellipses from the given parameters, but the following function should do the trick:

make_ellipse <- function(x0, y0, a, b, angle) {
  angle <- angle/360 * 2 * pi
  theta <- c(seq(0, 2 * pi, length.out = 360), 0)

  crds <- cbind(a * cos(theta) * cos(angle) - b * sin(theta) * sin(angle) + x0,
                a * cos(theta) * sin(angle) + b * sin(theta) * cos(angle) + y0)
  sf::st_polygon(list(crds))
}

This allows us to use your data frame to get the ellipses as follows:

a <- do.call(make_ellipse, data[1, 1:5])
b <- do.call(make_ellipse, data[2, 1:5])

And we can see that these are correct by doing:

plot(b)
plot(a, add = TRUE)

We can get the intersection as follows:

intersection <- sf::st_intersection(a, b)

Again, we can see this is correct by adding it to our plot:

plot(intersection, add = TRUE, col = 'red')

And since we are satisfied that the red intersection is indeed the shape whose area we are interested in, we can simply do:

sf::st_area(intersection)
#> [1] 105.3689

Answer 2

This is the same method as @Allan but using different packages: PlaneGeometry to deal with the two ellipses, and cgalPolygons (not on CRAN) to compute their intersection and the area. Good news: one finds almost the same result.

library(PlaneGeometry)
library(cgalPolygons)

ell1 <- Ellipse$new(center = c(0, 0), rmajor = 10, rminor = 5, alpha = 45)
ell2 <- Ellipse$new(center = c(0, 0), rmajor = 20, rminor = 3, alpha = 35)

ell1_path <- ell1$path(npoints = 200L)
ell2_path <- ell2$path(npoints = 200L)

ell1_polygon <- cgalPolygonWithHoles$new(ell1_path)
ell2_polygon <- cgalPolygonWithHoles$new(ell2_path)

# intersection
interList <- ell1_polygon$intersection(ell2_polygon)
inter <- interList[[1]]

# plot
opar <- par(mar = c(0, 0, 0, 0))
plot(
  NULL, xlim = c(-26, 26), ylim = c(-16, 16), asp = 1, 
  xlab = NA, ylab = NA, axes = FALSE
)
ell1_polygon$plot(list(lwd = 2), lwd = 2, density = 10, new = FALSE)
ell2_polygon$plot(list(lwd = 2), lwd = 2, density = 10, new = FALSE)
inter$plot(lwd = 4, col = "red", new = FALSE)
par(opar)

# area
inter$area()
# 105.357

---

Edit: ellipse partitioning

We can also split an ellipse and check if its area is close to the sum of the areas of its parts.

library(PlaneGeometry)
library(cgalPolygons)

# first ellipse
ell1 <- Ellipse$new(center = c(0, 0), rmajor = 10, rminor = 5, alpha = 45)
# second ellipse
ell2 <- Ellipse$new(center = c(0, 0), rmajor = 20, rminor = 3, alpha = 35)

# first ellipse as a path (2-columns matrix)
ell1_path <- ell1$path(npoints = 500L)
# second ellipse as a parh
ell2_path <- ell2$path(npoints = 500L)

# convert to polygons 
ell1_polygon <- cgalPolygonWithHoles$new(ell1_path)
ell2_polygon <- cgalPolygonWithHoles$new(ell2_path)

# intersection
ell12_polygon <- ell1_polygon$intersection(ell2_polygon)
# the intersection is a single polygon
ell12_polygon <- ell12_polygon[[1L]]

# plot the two ellipses and their intersection
opar <- par(mar = c(0, 0, 0, 0))
plot(
  NULL, xlim = c(-26, 26), ylim = c(-16, 16), asp = 1, 
  xlab = NA, ylab = NA, axes = FALSE
)
ell1_polygon$plot(list(lwd = 2), lwd = 2, new = FALSE)
ell2_polygon$plot(list(lwd = 2), lwd = 2, new = FALSE)
ell12_polygon$plot(lwd = 4, col = "red", new = FALSE)
par(opar)

# now subtract ell2 to ell1
subtraction_ell1_minus_ell2 <- ell1_polygon$subtract(ell2_polygon)
# this gives two polygons
plgn1 <- subtraction_ell1_minus_ell2[[1L]]
plgn2 <- subtraction_ell1_minus_ell2[[2L]]

# let's plot ell1 and its partition into three parts
opar <- par(mar = c(2, 2, 1, 1))
plot(
  NULL, xlim = c(-16, 16), ylim = c(-10, 10), asp = 1, 
  xlab = NA, ylab = NA, axes = FALSE
)
axis(1); axis(2)
# ellipse 1
ell1_polygon$plot(list(lwd = 2), lwd = 2, new = FALSE)
# ellipse 1 inter ellipse 2
ell12_polygon$plot(lwd = 2, col = "red", new = FALSE)
# the first piece of the subtraction
plgn1$plot(lwd = 2, col = "green", new = FALSE)
# the second piece of the subtraction
plgn2$plot(lwd = 2, col = "blue", new = FALSE)
par(opar)

# green area
garea <- plgn1$area()
# blue area
barea <- plgn2$area()
# red area
rarea <- ell12_polygon$area()

# sum of the three areas
garea + barea + rarea
# = 157.0755

# area of ellipse 1 (pi * major * minor)
pi * ell1$rmajor * ell1$rminor
# = 157.0796

The sum of the three areas is slightly less than the true area of the ellipse. That's normal because the paths we take approximate the true paths from below.

---

Edit: comparison with an "exact" result

There's an implementation on this Github repo. Here is an example treated with this implementation:

The resulting area is A12 and epsilon is the relative error.

I tried with my method. Here are the parameters:

ell1 <- Ellipse$new(
  center = c(14.527050, 8.343243), rmajor = 5.262357, rminor = 4.555325, alpha = -2.933767, degrees = FALSE
)
ell2 <- Ellipse$new(
  center = c(7.598473,  6.096601), rmajor = 6.720221, rminor = 5.767499, alpha = -1.624, degrees = FALSE
)

Then proceeding as above, taking 500 points to approximate the ellipses, I find 21.44292, very close to A12.