Biswajit Banerjee

### Plotting ellipsoids with R

Using the rgl graphics package.

- Introduction
- The
`rgl`

package - Drawing ellipsoids with
`cda`

- Euler angles and direction cosines
- Plotting ellipsoids with
`rgl.ellipsoid`

- Plotting ellipsoids with
`my.ellipsoid`

- Plotting all the ellipsoids with
`my.ellipsoid`

- Remarks

#### 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

and then load the package in our R script using

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:

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,

The function to draw N ellipsoids is called using

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.

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:

Here `R`

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

##### The image produced by `cda`

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

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

The Euler angles can be converted into direction cosines using where

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

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

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.

The R function used to produce these images is given below

and the function is called using

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:

and used this function in my ellipsoid drawing function:

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.

#### Plotting all the ellipsoids with `my.ellipsoid`

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.

This function was called as shown below.

#### 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.

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

📅 06.11.2017 📁 R