# Airy Area: approximating surface area of a cell from a 3D point set

**Rstats – quantixed**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.

**leave a comment**for the author, please follow the link and comment on their blog:

**Rstats – quantixed**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.