Biswajit Banerjee

Plotting ellipsoids with R

Using the rgl graphics package.

Introduction

Recently I ran into the problem of having to visualize the locations and orientations of around 700 ellipsoidal particles. The input data were in a CSV text file and I decided to use R to process and plot the data.

R has powerful file parsing packages that makes reading the data trivial. However, I wasn’t sure of the 3D graphics capabilities of R.

In this article, I will talk about my experience plotting 3D graphics in R.

The rgl package

R has a useful OpenGL-based graphics package called rgl. To install the package we use

install.packages('rgl')

and then load the package in our R script using

require('rgl')

To check that the rgl package is producing the right output, we can draw axis arrows on the render window. Assume that the axes are aligned along , , and and have unit length. The following script produces the image below:

open3d()
arrow3d(c(0, 0, 0), c(1, 0, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="red")
arrow3d(c(0, 0, 0), c(0, 1, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="green")
arrow3d(c(0, 0, 0), c(0, 0, 1), s = 1/7, width = 1/5, n = 10, type="rotation", col="blue")
axes3d()
rgl.snapshot("rgl_axis.png")
RGL axis arrows
Figure 1. Axes generated with rgl. The image has been trimmed with ImageMagick using the command
"mogrify -trim rgl_axis.png"

The arrow3d function takes the start and end points of the arrow, a parameter s that controls the size of the arrowhead, the width parameter that controls the arrow thickness, n to control the number of facets in the cylinder and head of the arrow, a type parameter that determines the shape of the arrow, and material parameters such as color.

The axes3d function generates the bounding box with tick marks and text. To save the image, we use the rgl.snapshot command.

Drawing ellipsoids with cda

While ellipsoids can be drawn directly using the ellipse3d function in rgl, there is a lot of extra work done by that function that is not needed for our purposes. I decided to try the simpler version called rgl.ellipsoids provided by the cda package in R. As before, we install and load the package using,

install.packages('cda')
require('cda')

The function to draw N ellipsoids is called using

rgl.ellipsoids(pos_mat, size_mat, euler_angles, col="gold")

where pos_mat is a 3 x N matrix of ellipsoid centers, size_mat is a 3 x N matrix of ellipsoid semi-axes, and euler_angles is a 3 x N matrix of Euler angles in the sequence . Our input data on the orientations of the ellipsoids was in the form of angles with respect to the reference coordinate axes. These has to be converted into direction cosines and then into Euler angles.

Euler angles using orientlib

To do the conversion, I first tried to search for a package that can do the conversion. I found that the orientlib package appeared to be capable of giving me what I needed and so I installed and loaded that.

install.packages('orientlib')
require('orientlib')

To get the Euler angles from orientlib, we have to reshape the direction cosines from out input data into a 3 x 3 matrix whose columns are the direction cosines of the three axis directions of the ellipsoid. We use the ZXZ convention, and use the following commands to generate the Euler angles:

Rmat <- rotmatrix(t(R))
angles <- eulerzxz(Rmat)
phi_a = angles@x[1,1]
theta_a = angles@x[1,2]
psi_a = angles@x[1,3]

Here R is the rotation matrix whose rows are the direction cosines of the three axes of an ellipsoid.

The image produced by cda
RGL axis arrows
Figure 2. Ellipsoids generated with cda. The orientations of the ellipsoids are not correct because some of them appear to intersect.

The cda manual does not talk about the convention used to compute the rotation matrix from the input Euler angles. If we assume the ZXZ convention, we get erroneous results because ellipsoid appear to intersect as can be seen in Figure 2. We also find that the axis arrows are removed during the process of drawing the ellipsoids.

As we were unsure whether the problem was with the Euler angles computed by orientlib or with the rotation matrices computed by cda, we decided to implement the Euler angle calculation ourselves.

Euler angles and direction cosines

The calculation of Euler angles is straightforward if we assume the ZXZ convention. In particular, we adopt the convention that the direction cosine matrix, also called the rotation matrix (), is given by

$$ \mathbf{R} = \begin{bmatrix} \mathbf{E}_x \cdot \mathbf{e}_a & \mathbf{E}_y \cdot \mathbf{e}_a & \mathbf{E}_z \cdot \mathbf{e}_a \\ \mathbf{E}_x \cdot \mathbf{e}_b & \mathbf{E}_y \cdot \mathbf{e}_b & \mathbf{E}_z \cdot \mathbf{e}_b \\ \mathbf{E}_x \cdot \mathbf{e}_c & \mathbf{E}_y \cdot \mathbf{e}_c & \mathbf{E}_z \cdot \mathbf{e}_c \end{bmatrix} = \begin{bmatrix} a_x & a_y & a_z \\ b_x & b_y & b_z \\ c_x & c_y & c_z \end{bmatrix} $$

where are unit vectors along the reference coordinate axes and are unit vectors along the axes of the ellipsoid.

Then, the Euler angles are given by

$$ \begin{align} \theta & = \cos^{-1}(R_{33}) \\ \phi & = \begin{cases} \tan^{-1}(R_{21}/R_{11}) & \quad \text{for} \quad \theta \sim 0 \\ \tan^{-1}(-R_{13}/R_{23}) & \quad \text{otherwise} \end{cases}\\ \psi & = \begin{cases} 0 & \quad \text{for} \quad \theta \sim 0 \\ \tan^{-1}(R_{31}/R_{32}) & \quad \text{otherwise} \end{cases} \end{align} $$

The Euler angles can be converted into direction cosines using where

$$ \mathbf{D} = \begin{bmatrix} \cos\phi & -\sin\phi & 0 \\ \sin\phi & \cos\phi & 0 \\ 0 & 0 & 1 \end{bmatrix} ~,~~ \mathbf{C} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \end{bmatrix} ~,~~ \mathbf{B} = \begin{bmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

We can easily check that these transformations produce orthogonal rotation matrices.

isRotationMatrix <- function(R) {
  Rt <- t(R)
  shouldBeIdentity <- Rt %*% R
  I <- diag(3)
  n <- norm(I - shouldBeIdentity, "2")
  return(n < 1.0e-6)
}

The R function that we use to do the Euler angle computations is given below.

rotationMatrixToEulerAngles <- function(R) {
  # Compute Euler angles
  costheta = R[3,3]
  sinthetaSq = 1 - costheta^2
  if (sinthetaSq < 1.0e-20) {
    theta = acos(costheta)
    phi = atan2(R[2,1], R[1,1])
    psi = 0
  } else {
    theta = acos(costheta)
    phi = atan2(R[1,3], -R[2,3])
    psi = atan2(R[3,1],  R[3,2])
  }
  # Check rotation matrix
  cp = cos(phi)
  sp = sin(phi)
  ct = cos(theta)
  st = sin(theta)
  cs = cos(psi)
  ss = sin(psi)
  D = matrix(c(cp, -sp, 0, sp, cp, 0, 0, 0, 1), nrow = 3, ncol = 3, byrow = TRUE)
  C = matrix(c(1, 0, 0, 0 , ct, -st, 0, st, ct), nrow = 3, ncol = 3, byrow = TRUE)
  B = matrix(c(cs, -ss, 0, ss, cs, 0, 0, 0, 1), nrow = 3, ncol = 3, byrow = TRUE)
  DCB = D %*% C %*% B
  if (!isRotationMatrix(R) || !isRotationMatrix(DCB)) {
    diff = DCB - R
    print(diff)
  }
  # Return Euler angles
  return(c(phi, theta, psi))
}

Euler angles computed using the above function are identical to those produced by orientlib.

Plotting ellipsoids with rgl.ellipsoid

We now tried to plot individual ellipsoids using the rgl.ellipsoid function provided by cda. We noticed that the orientations of the largest axis of some ellipsoids and the shapes produced by cda were not aligned.

RGL axis arrows RGL axis arrows
Figure 3. Ellipsoids 2 an 3 generated with cda. The orientation of ellipsoid 2 appears correct but that of ellipsoid 3 is not.

The R function used to produce these images is given below

plotCDAEllipsoid <- function(index, pos, size, euler, dircos) {
  va = dircos[1:3]
  vb = dircos[4:6]
  vc = dircos[7:9]
  ppa = pos + size[1]*va
  ppb = pos + size[2]*vb
  ppc = pos + size[3]*vc
  arrow3d(pos, ppa, type="flat", col="red")
  arrow3d(pos, ppb, type="flat", col="green")
  arrow3d(pos, ppc, type="flat", col="blue")
  ell1 = rgl.ellipsoid(x = pos[1], y = pos[2], z = pos[3],
                       a = size[1], b = size[2], c = size[3],
                       phi = euler[1], theta = euler[2], psi = euler[3],
                       subdivide = 3, smooth = TRUE, color = "gold", alpha=0.5)
  plot3d(ell1, add = TRUE)
  text3d(x = pos[1], y = pos[2], z = pos[3], as.character(index))
}

and the function is called using

open3d()
for (i in c(2,3)) {
  plotCDAEllipsoid(i, pos_mat[,i], size_mat[,i], euler_mat[,i], cos_angles[i,])
}
arrow3d(c(0, 0, 0), c(0.0006, 0, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="red")
arrow3d(c(0, 0, 0), c(0, 0.0006, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="green")
arrow3d(c(0, 0, 0), c(0, 0, 0.0006), s = 1/7, width = 1/5, n = 10, type="rotation", col="blue")
axes3d()
view3d(fov=0)

Notice that an orthogonal projection is used to view the result (fov = 0) so that distortions are avoided.

Plotting ellipsoids with my.ellipsoid

After some exploration, it became clear that my definitions of did not match those used by cda. So I decided to rewrite the cda functions to use the rotation matrix directly:

my.ellipsoid <- function (pos, size, euler, subdivide = 3, smooth = TRUE, ...)
{
  x = pos[1]
  y = pos[2]
  z = pos[3]
  a = size[1]
  b = size[2]
  c = size[3]
  phi = euler[1]
  theta = euler[2]
  psi = euler[3]
  cp = cos(phi)
  sp = sin(phi)
  ct = cos(theta)
  st = sin(theta)
  cs = cos(psi)
  ss = sin(psi)
  D = matrix(c(cp, -sp, 0, sp, cp, 0, 0, 0, 1), nrow = 3, ncol = 3, byrow = TRUE)
  C = matrix(c(1, 0, 0, 0 , ct, -st, 0, st, ct), nrow = 3, ncol = 3, byrow = TRUE)
  B = matrix(c(cs, -ss, 0, ss, cs, 0, 0, 0, 1), nrow = 3, ncol = 3, byrow = TRUE)
  rotMat = D %*% C %*% B

  sphere <- rgl::subdivision3d(rgl::cube3d(...), subdivide)
  class(sphere) <- c("mesh3d","shape3d")

  norm <- sqrt(sphere$vb[1, ]^2 +
               sphere$vb[2, ]^2 +
               sphere$vb[3, ]^2 )
  for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm
  sphere$vb[4, ] <- 1
  sphere$normals <- sphere$vb
  result <- rgl::scale3d(sphere, a,b,c)
  result <- rgl::rotate3d(result,matrix=rotMat)
  result <- rgl::translate3d(result, x,y,z)
  invisible(result)
}

and used this function in my ellipsoid drawing function:

plotMyEllipsoid <- function(index, pos, size, euler, dircos) {
  va = dircos[1:3]
  vb = dircos[4:6]
  vc = dircos[7:9]
  ppa = pos + size[1]*va
  ppb = pos + size[2]*vb
  ppc = pos + size[3]*vc
  arrow3d(pos, ppa, type="flat", col="red")
  arrow3d(pos, ppb, type="flat", col="green")
  arrow3d(pos, ppc, type="flat", col="blue")

  ell2 = my.ellipsoid(pos, size, euler, subdivide = 3, smooth = TRUE, color = "gold", alpha=0.5)
  plot3d(ell2, add = TRUE)
  text3d(x = pos[1], y = pos[2], z = pos[3], as.character(index))
}
RGL my ellipsoids 2 and 3
Figure 4. Ellipsoids 2 an 3 generated with my ellipsoid plotter. The orientations of both ellipsoids are now correct.

As with the CDA ellipsoid plotter, my ellipsoid plotter is called using similar code. The resulting plot produces the correct orientations, as can be seen in Figure 4. If particular, notice that the red arrows inside the ellipsoids, which correspond to the longest axis of each ellipsoid, are align correctly with the plotted ellipsoids as do the other two axes.

open3d()
for (i in c(2,3)) {
  plotMyEllipsoid(i, pos_mat[,i], size_mat[,i], euler_mat[,i], cos_angles[i,])
}
arrow3d(c(0, 0, 0), c(0.0006, 0, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="red")
arrow3d(c(0, 0, 0), c(0, 0.0006, 0), s = 1/7, width = 1/5, n = 10, type="rotation", col="green")
arrow3d(c(0, 0, 0), c(0, 0, 0.0006), s = 1/7, width = 1/5, n = 10, type="rotation", col="blue")
axes3d()
view3d(fov=0)

Plotting all the ellipsoids with my.ellipsoid

RGL my ellipsoids 2 and 3
Figure 5. Ellipsoids generated with my ellipsoid plotter. Note that there are no overlapping ellipsoids.

Now that we have figured out how to plot the ellipsoid in the correct orientation, the next step is to plot the full set more efficiently using my version of rgl.ellipsoids called my.ellipsoids. Note that we are adding the axis arrows for the plot in this function and creating a list of shapes to pass into shapelist3d.

A plot of the ellipsoids produced by this approach is shown in Figure 5. Note that now there are no longer any overlaps.

my.ellipsoids <- function(positions, sizes, angles, colour = "red", ...){
  N <- NCOL(positions)
  colours <- rep(colour, length.out=N)
  ll <- lapply(seq(1,N), function(ii)
    my.ellipsoid(positions[,ii], sizes[,ii], angles[,ii],, col = colours[ii], ...))

  a1 <- arrow3d(c(0.002, 0, 0), c(0.0025, 0, 0), s = 1/3, width = 1/3, n = 10, type="rotation", col = "red", plot = FALSE)
  a2 <- arrow3d(c(0.002, 0, 0), c(0.002, 0.0005, 0), s = 1/3, width = 1/3, n = 10, type="rotation", plot = FALSE, col="green")
  a3 <- arrow3d(c(0.002, 0, 0), c(0.002, 0, 0.0005), s = 1/3, width = 1/3, n = 10, type="rotation", plot = FALSE, col="blue")

  rgl::shapelist3d(c(ll, list(a1, a2, a3)))
}

This function was called as shown below.

open3d()
my.ellipsoids(pos_mat, size_mat, euler_angles, col="gold", add=TRUE)
axes3d()
view_matrix = matrix(c(-0.21751294, -0.53017938,  0.81950366,  0.00000000,  
                        0.97544736, -0.08853958,  0.20162290,  0.00000000,
                       -0.03433822,  0.84324223,  0.53642339,  0.00000000,
                        0.00000000,  0.00000000,  0.00000000,  1.00000000),
                     nrow = 4, ncol = 4)
view3d(userMatrix = view_matrix)
rgl.snapshot("rgl_my_ellipsoids.png")

Remarks

The above exercise shows that once can get reasonably quick 3D graphics plots of ellipsoids using R. For completeness, I have listed below how the data were read into R.

#------------------------------------------------------------------
# 1) Read the particle data
#------------------------------------------------------------------
partDistCSV <- "input_particle_file_orig"

# Read the first line (contains the number of particles)
df_numPart <- read.csv(partDistCSV, header = FALSE, nrows = 1, sep="",
                   blank.lines.skip = TRUE)
numPart <- df_numPart$V1

# Read the particle data into a data frame
# Data should contain the following header:
# > names(df)
# [1] "id"         "type"       "radius_a"   "radius_b"   "radius_c"  
# [6] "position_x" "position_y" "position_z" "axle_a_x"   "axle_a_y"  
#[11] "axle_a_z"   "axle_b_x"   "axle_b_y"   "axle_b_z"   "axle_c_x"  
#[16] "axle_c_y"   "axle_c_z"   "velocity_x" "velocity_y" "velocity_z"
#[21] "omga_x"     "omga_y"     "omga_z"     "force_x"    "force_y"   
#[26] "force_z"    "moment_x"   "moment_y"   "moment_z"  
df_part <- read.csv(partDistCSV, header = TRUE, skip = 2, sep="",
                   blank.lines.skip = TRUE, nrows = numPart)

#------------------------------------------------------------------
# 2) Create positions and radii
#------------------------------------------------------------------
positions <- data.frame(x = df_part$position_x, y = df_part$position_y, z = df_part$position_z)
sizes <- data.frame(x = df_part$radius_a, y = df_part$radius_b, z = df_part$radius_c)
pos_mat <- data.matrix(t(positions))
size_mat <- data.matrix(t(sizes))

#------------------------------------------------------------------
# 3) Create Euler angles
#------------------------------------------------------------------
angles <- data.frame(xa = df_part$axle_a_x, ya = df_part$axle_a_y, za = df_part$axle_a_z,
                     xb = df_part$axle_b_x, yb = df_part$axle_b_y, zb = df_part$axle_b_z,
                     xc = df_part$axle_c_x, yc = df_part$axle_c_y, zc = df_part$axle_c_z)
cos_angles <- cos(angles)
cos_mat <- data.matrix(cos_angles)
euler_angles <- apply(cos_angles, 1,
                 function(x) {
                   R <- matrix(x, nrow = 3, ncol=3, byrow=TRUE)
                   angles = rotationMatrixToEulerAngles(R)
                   return(angles)
                 })
euler_mat <- data.matrix(euler_angles)

If you have questions/comments/corrections, please contact banerjee at parresianz dot com dot zen (without the dot zen).

 