In the spirit of “if it took you a while to find out how to do something, write about it”, I will detail a method to approximate the surface area of a 3D shape.

Our application here was finding the surface area of a cell but it can be used on any shape.

We start with a 3D point set, specified by points of interest in a single cell imaged by 3D confocal microscopy. The details of the points aren’t too important, but they are generally on the cell surface with some internal points.

This is an ugly plot, don’t worry, it will get better!

To find the surface of this point set we can use 3D alpha shapes. Briefly, this is a method to trim away at a 3D convex hull to more finely contour the surface of an object. Fortunately, there is a nice R package called alphashape3d on CRAN to do the calculation.

library(alphashape3d)
# Create a data frame
df <- read.csv("Data/measure3d2.csv")
# the x y z coords are in columns 8 to 10 of the file, stored in real units (um)
mata <- as.matrix(df[,8:10])
# 3d coords can be fed into the function ashape3d()
# the second argument specifies alpha
alfa <- ashape3d(mata, 5)
plot(alfa)

the plot() function here displays the shape using {rgl} and XQuartz

So we have our 3D shape. There is a function in {alphashape3d} to calculate volume:

volume_ashape3d(alfa)
>2151.8822004701

The question is: how can we calculate the surface area?

I couldn’t see a way to do this using an in-built function, so I rolled my sleeves up.

The solution

We simply need to sum all the triangles represented in the plot above to get an approximation of the surface.

If we take a look at alfa, we can see it is an S3 object termed ashape3d. View(alfa)

The matrix triang contains the triangle information, while x contains the 3D point set that we used to generate the ashape3d. From the docs we can see that the ninth column of triang tells us whether that triangle is in the alpha shape or not (if it has a value or 2 or 3). We also know that the columns termed tr1, tr2 and tr3 specify the indices of points (of x) that make up the triangle.

Putting that information together, we can do something like:

selection <- alfa$triang[, 9] >= 2
triangles <- alfa$triang[selection, c("tr1", "tr2", "tr3")]
triangles <- as.data.frame(triangles)
triangles$area <- calculate_area(alfa$x, triangles)

this extracts the triangles of interest and collects the three points that comprise them. Then we send this, together with the point set, to a function in order to do the calculation.

This function, calculate_area() uses a bunch of subfunctions to get the area of the triangles of interest. At the heart of it is Heron’s formula.

# function to calculate area of a triangle in 3D space
heron <- function(a,b,c){
  s <- (a + b + c) / 2
  area <- sqrt(s*(s-a) * (s-b)*(s-c))
  return(area)
}

# utility to get the distance between two 3D coords
distance3d <- function(x1,y1,z1,x2,y2,z2){
  a <- (x1-x2)^2+(y1-y2)^2 + (z1-z2)^2
  d <- sqrt(a)
  return(d)
}

# wrapper to calculate the area of triangle from 2 x 3d coord sets
areatriangle3d <- function(x1,y1,z1,x2,y2,z2,x3,y3,z3){
  a <- distance3d(x1,y1,z1,x2,y2,z2)
  b <- distance3d(x2,y2,z2,x3,y3,z3)
  c <- distance3d(x3,y3,z3,x1,y1,z1)
  A <- heron(a,b,c)
  # print(paste("area of triangle is",A))
  return(A)
}

# this function calculates the area of each triangle
calculate_area <- function(x, triangles) {
  # x is a matrix of 3D coordinates
  # triangles is a data frame of triangle indices
  # calculate the area of each triangle
  
  a <- areatriangle3d(x[triangles$tr1, 1], x[triangles$tr1, 2], x[triangles$tr1, 3],
                      x[triangles$tr2, 1], x[triangles$tr2, 2], x[triangles$tr2, 3],
                      x[triangles$tr3, 1], x[triangles$tr3, 2], x[triangles$tr3, 3])
  return(a)
}

Now triangles has a column containing the area of each triangle of interest, and we can simply sum them to approximate the cell surface.

>sum(triangles$area)
1107.63245400761

Well, there’s a solution, or perhaps the solution, to the problem. But how do we know that it’s correct?

Verifying the method

Let’s use a simple shape to check whether the method works.

# construct a 9 x 3 matrix of points on the surface of a cube
cube <- matrix(c(-1,1,-1,1,-1,1,-1,1,
         1,1,-1,-1,1,1,-1,-1,
         1,1,1,1,-1,-1,-1,-1), ncol = 3)
cube_a <- ashape3d(cube, 10)
plot(cube_a)

And this gives:

Which is a cube with edges of 2 units, centred on the origin.

We can use a wrapper to do the calculations for us:

# this functions does all the calculations and prints the results
calculations_3d <- function(mat) {
  
  alphamat <- ashape3d(mat, 5)
  
  # volume of alpha shape
  cat(paste0("Volume: ",volume_ashape3d(alphamat),"\n"))
  # surface area of alpha shape
  selection <- alphamat$triang[, 9] >= 2
  triangles <- alphamat$triang[selection, c("tr1", "tr2", "tr3")]
  triangles <- as.data.frame(triangles)
  triangles$area <- calculate_area(alphamat$x, triangles)
  # sum the area of all triangles
  cat(paste0("Surface area: ",sum(triangles$area),"\n"))
  # number of spots that went in
  cat(paste0("Number of spots: ",nrow(mat),"\n"))
  # density of spots is the number of spots divided by the surface area
  cat(paste0("Density of spots: ",nrow(mat) / sum(triangles$area),"\n"))
  # how many spots are in the alpha shape? the number of unique triangle indices
  cat(paste0("Number of spots in alpha shape: ",length(unique(c(triangles$tr1, triangles$tr2, triangles$tr3))),"\n"))
  # density of the spots in alpha shape
  cat(paste0("Density of spots in alpha shape: ",length(unique(c(triangles$tr1, triangles$tr2, triangles$tr3))) / sum(triangles$area),"\n"))
}

So, let’s try it:

> calculations_3d(cube)
Volume: 8
Surface area: 24
Number of spots: 8
Density of spots: 0.333333333333333
Number of spots in alpha shape: 8
Density of spots in alpha shape: 0.333333333333333

The surface area is indeed 24 units, as expected.

Let’s try a different object: a sphere.

We first need a 3D point set. I found this handy function online. We can make a sphere with a radius of 1, and either have points throughout the sphere or just on the surface.

rsphere <- function(n, r = 1.0, surface_only = FALSE) {
  phi       <- runif(n, 0.0, 2.0 * pi)
  cos_theta <- runif(n, -1.0, 1.0)
  sin_theta <- sqrt((1.0-cos_theta)*(1.0+cos_theta))
  radius <- r
  if (surface_only == FALSE) {
    radius <- r * runif(n, 0.0, 1.0)^(1.0/3.0)
  }
  
  x <- radius * sin_theta * cos(phi)
  y <- radius * sin_theta * sin(phi)
  z <- radius * cos_theta
  
  cbind(x, y, z)
}

# when we test the function, we get the following results
set.seed(123)
sphere_points <- rsphere(500, surface_only = T)
calculations_3d(sphere_points)

# plot the sphere
plot3d(sphere_points, col = "red", size = 0.5)
sphere_alpha <- ashape3d(sphere_points, 10)
plot(sphere_alpha)

The results:

Volume: 4.09187028812644
Surface area: 12.4200231832229
Number of spots: 500
Density of spots: 40.2575738083488
Number of spots in alpha shape: 500
Density of spots in alpha shape: 40.2575738083488

The volume is approximately 4/3 * pi * r^3 = 4.19 and the surface area is approximately 4 * p i* r^2 = 12.57

If we increase the number of points we get closer to these numbers, although the method using triangles will always be less than the true value.

> set.seed(123)
> sphere_points <- rsphere(5000, surface_only = T)
> calculations_3d(sphere_points)
Volume: 4.17880848901202
Surface area: 12.5513891616036
Number of spots: 5000
Density of spots: 398.362279714479
Number of spots in alpha shape: 5000
Density of spots in alpha shape: 398.362279714479

Finally, these checks have only used the situation of points that are on the surface of the object, let’s make some internal points:

> set.seed(123)
> sphere_points <- rsphere(10000, surface_only = F)
> calculations_3d(sphere_points)
Volume: 4.00201444091866
Surface area: 12.2429156824266
Number of spots: 10000
Density of spots: 816.798894919607
Number of spots in alpha shape: 471
Density of spots in alpha shape: 38.4712279507135

In this example, we generated 10000 points anywhere in the sphere, but only 471 were on the alpha shape, rather than being internal. Our approximation is similar to the 500 point (all surface) result above.

Conclusion

Even though there was no in-built function to approximate the surface area of a 3D alpha shape. It is possible to do this direct from the ashape3d object.

—

The post title comes from “Airy Area” by Ozric Tentacles from the “Vitamin Enhanced” album.