Zoom and rotate the ammonite (from mobile, only zooming is allowed)

   

Introduction

Among all life forms, both living and extinct, few embody geometry as elegantly as ammonites. Their beauty crosses the boundaries of palaeontology and is inherently satisfying to the eye of both the fossil expert and the layman. I am neither a layman nor a fossil expert (I am somewhere in between), yet I go fossil-hunting and the best among my finds is frustration. I long the retrieval of intact ammonite shells but my pockets are room for only fragments. My skills, though, are not those of a fossil-hunter, so I decided to become a fossil manufacturer, at least from the comfort of my laptop. By using simple DIY geometry and some tricks, here I will explain how to produce your very own collection of virtual ammonites.

(If you are itching to generate ammonites in 3D and want to skip the how to part, go to the bottom of this document.)

 

Step 1: from circle to spiral

To build our ammonite, first we start with a circle. Let’s work in 2D (we will move to the third dimension later on). Building a circle in R is easy peasy. Let’s define a number of points along the perimeter of the circle (n=50), then we need a radius (r=1) and sequential angles (equally spaced) for positioning the points around the circle. The angles go from 0 to 2*pi (in degrees, that would be from 0 to 360). Finally, we need to specify the x and y of the center of the circle (ctr):

n<-50
r<-1
ctr<-c(0,0)

th<-seq(0,2*pi,length=n)

x<-ctr[1]+r*cos(th)
y<-ctr[2]+r*sin(th)

plot(x,y,pch=19,cex=0.8,asp=1)

The x and y coordinates of the circle are then calculated using ctr[1]+r*cos(th) and ctr[2]+r*sin(th), respectively (for simplicity, we set the centre of the circle at 0, but we could use any xy point).

An ammonite shell is not a circle, though, but a spiral. A spiral is nothing else than a circle whose ends don’t like each other and try keeping apart while rotating over and over. To make this happen in R, we need to use increasing radii for each consecutive point. In particular, we want to produce a logarithmic spiral, with radii increasing exponentially from the center to the outer points. The radii have to move around the center more than 2*pi, because the spiral goes around the centre multiple times, which we set as turns (we also set some new larger n):

n<-300
turns<-20
th<-seq(0,turns*pi,length=n)

Now that we have the angles, we set the radii of each point lying at said angles. To do so, we use two parameters, a and b, that modify the morphology of the spiral through size and rate of radii increase, respectively. When the parameter b increases, the spiral move outward faster, pretty much like in the shell of the living Nautilus. The formula for the radii of a spiral is

\[ r = ae^{b\theta} \] where \(\theta\) is the angle (th in our code):

a<-1
b<-0.13

r<-a*exp(b*th)

x<-ctr[1]+r*cos(th)
y<-ctr[2]+r*sin(th)

The code for generating the points along the spiral is the same we used for the circle (you can change the values of a and b but be careful with the latter - keep b between 0.1 and 0.3). Now you can plot the spiral:

plot(x,y,cex=0.5,asp=1,type="l",pch=19)

 

Step 2: make it 3D

It’s now time to move our spiral into the third dimension. What we need is a zero, nothing more. We have already computed x and y, now we need z, and we set it as 0:

z<-0
spi<-cbind(x,y,z)

Let’s put x, y and z into one single matrix (spi), which includes the coordinates of our 3D spiral. We can plot it using the R package rgl. You can use the graph below to move it and look at the skeleton of our ammonite:

library(rgl)
lines3d(spi,lwd=2)

 

Step 3: give it some depth

We have produced a logarithmic spiral whose morphology depends on the parameters a, b and turns. How do we go from this template to the actual ammonite, with its three-dimensional attributes?

The spiral is composed of 3D points and at each of these points we give some depth in the form of a section. In other words, we create the shape of the ammonite in the third dimension and, for simplicity, we make the section circular. We define nsection as the number of 3D points of each circular section, and generate the angles for the nsection points of each section:

nsection<-30

th<-seq(0,2*pi,length=nsection+1)
th<-th[-length(th)]

Each section will have to decrease in size from the outermost point to the centre of the spiral, and that needs another parameter defining the thickness of the spiral at each point. We call this parameter thick and we use it to reduce the radius of each section proportionally to its position along the spiral (the result is the object rsection)

thick<-2.5

rsection<-r/thick

Let’s see how to position the sections along the spirals by looking at one example (then we will proceed to compute all the sections). Above, we have seen how to compute the coordinates of a circle, and we do it again here for generating the circular section. Remember that our spiral lies on the xy plane, with z=0. For our circular section to be perpendicular to the plane of the spiral, we need to set y=0. As the center of the section, let’s use the outermost point (the last of the object spi); also, the radius we use is not r but rsection:

x<-spi[n,1]+rsection[n]*cos(th)
y<-spi[n,2]
z<-spi[n,3]+rsection[n]*sin(th)

sect<-cbind(x,y,z)

Let’s plot the spiral and the section to see what we have done so far:

lines3d(spi,lwd=2)
lines3d(sect,col="red",radius=2)

The plane of the section is perpendicular to the plane of the spiral but we need to rotate it to give it the correct orientation. The orientation is found by rotating the section around its z axis (in other words, it rotates but the plane of the section stays perpendicular to the plane of the spiral).

Bear with me here. For the rotation to occur, we need to know what axis it will rotate around (we already know that is the z axis) and we need to know of what angle to rotate it. This angle is the angle between the normal to the section (the direction that is perpendicular to the section plane) and the direction from one spiral point to the next. If this does not make much sense to you, look at the following plot. The red and blue circles are the original and rotated sections, respectively; the red arrow is the normal vector to the plane of the original section and the blue arrow is the vector direction from one point of the spiral to the next. What we want is to find the angle between the two arrows and use it to rotate the red section into the blue one.

The calculation of the angle takes a few steps:

#Calculate vector direction from the point on the spiral and the first point of the circuar section. The difference between two points is a direction! (Remember that, as an example, we are considering only one section here, the outermost one, whose spiral point we can collect using `n`)
a<-spi[n,]-sect[1,]

#Do the same but using another direction (let's take the one that lies at about 90 degrees from the one previously calculated)
b<-spi[n,]-sect[round(nsection/4),]

#The code below computes the cross product of the two vectors previously found: the result is a new vector that is perpendicular to both of them. The resulting vector is the normal to the plane of the section (perpendicular to such plane)
m<-outer(a,b,"*")
nor<-c(m[2,3]-m[3,2],
       m[3,1]-m[1,3],
       m[1,2]-m[2,1])

#We need the normal to the section because we want to calculate the angle existing between its normal and the direction existing between the point of the spiral (which is the centre of the section) and the next point on the spiral. Here, we compute this vector direction
dif<-spi[n,]-spi[n-1,]

#Now that we obtained the normal to the section plane and the vector direction from one point of the spiral to the next, we can compute the angle between the two vectors. This is the angle that we will use for rotating the section around the z axis. We add an 'if' statement because the vector direction from consecutive points of the spiral changes depending on our position around the spiral (based on the sign of the coordinates in 'dif')
ang<-acos(nor %*% dif/(sqrt(nor %*% nor) * sqrt(dif %*% dif)))
if(dif[1]<0){ang<- -ang}

Now that we have the angle, we can rotate the section around the z axis of an angle equal to ang. To do so, we use the function rotate3d from the package rgl. Since the rotation is referenced to the origin of the axes (point x=0, y=0, z=0), the section has to be translated back onto the spiral. We do so using the rgl function translate3d, but first we calculate the direction of the translation (tr):

rot<-rotate3d(sect,-c(ang),0,0,1)

tr<-spi[n,]-apply(rot,2,mean)
rot<-translate3d(rot,tr[1],tr[2],tr[3])

All these steps get your section rotated so that its normal vector is in the same direction to the direction between its center and the next point of the spiral.

 

Step 4: build all sections

We have seen the operations we need to perform for each section to reach the correct orientation. Let’s apply them to all section at each point of the spiral. This can be easily achieved using a for loop. The points of each section correctly oriented are stored in the object pts:

pts<-NULL
for(i in nrow(spi):2){
  x<-spi[i,1]+rsection[i]*cos(th)
  y<-spi[i,2]
  z<-spi[i,3]+rsection[i]*sin(th)
  sect<-cbind(x,y,z)
  
  a<-spi[i,]-sect[1,]
  b<-spi[i,]-sect[round(nsection/4),]
  m<-outer(a,b,"*")
  nor<-c(m[2,3]-m[3,2],
         m[3,1]-m[1,3],
         m[1,2]-m[2,1])
  
  dif<-spi[i,]-spi[i-1,]
  
  ang<-acos(nor %*% dif/(sqrt(nor %*% nor) * sqrt(dif %*% dif)))
  if(dif[1]<0){ang<- -ang}
  
  rot<-rotate3d(sect,-c(ang),0,0,1)
  tr<-spi[i,]-apply(rot,2,mean)
  rot<-translate3d(rot,tr[1],tr[2],tr[3])
  
  pts<-rbind(pts,rot)
}

Let’s visualise all sections together, which will give us an idea of the final ammonite:

lines3d(spi,col="red",lwd=2)
points3d(pts,size=0.5)

 

Step 5: triangulation

Our ammonite is now a point cloud, which we can use to build a 3D mesh. A mesh is more than just its points: it includes the information of how the points are connected to each other. We will build a triangular mesh, meaning that each point will be connected to others to form triangles, which are the faces of the mesh. Each triangle is expressed as the index of the points that are connected to form a triangle (the index is the position of the points in the object pts).

The triangulation can be performed using several algorithms (e.g. Delaunay algorithms). Instead, we will produce a bespoke, rather flawless triangulation. We can do that because each section is formed by the same number of points (here, we used nsection=30) laid out sequentially along the section, and each section is oriented in the same way with respect to the spiral (i.e. the point 5 of one section will always face the point 5 of the adjacent section). Based on these premises, the point 1 of the first section will be facing the point 1+nsection (point 31 of the object pts) of the following section, and the closest point with which a triangle can be formed is, for example, 32 (because it is adjacent to the point 31). We can generalise this rule and find all the triangles using the following code:

#Generate the indices of the positions for each point in the object 'pts'
pos<-1:nrow(pts)

#Associate at each point index, the facing point in the following section ('nsection+pos') and the point adjacent to that one ('nsection+pos+1'). We also make sure that our operations do not introduce indices higher than the maximum number of points
tr1<-cbind(pos,nsection+pos,nsection+pos+1)
  tr1<-ifelse(tr1>max(pos),NA,tr1)
  tr1<-na.omit(tr1)
  
#We do the same for the other haf of triangles, which are found by associating to each point its neighbour in the section ('pos+1') and the facing point in the following section ('nsection+pos+1)
tr2<-cbind(pos,pos+1,nsection+pos+1)
  tr2<-ifelse(tr2>max(pos),NA,tr2)
  tr2<-na.omit(tr2)

#This is the first triangle, which I could not include in the rule above (not an elegant solution but it works)
tr0<-c(1,nsection,nsection+1)
  
#We bind all trinagle indices together in one matrix
tri<-rbind(tr0,tr1,tr2)

To be able to use these information as a 3D mesh, we need to put them together in a suitable format. In R, a common format for 3D meshes is the class mesh3d, used also by the package rgl. This format is a list containing the points and the triangles (besides other attributes that we will ignore here), but we need to make some alterations to them to make it work:

mesh<-mesh3d(vertices=t(cbind(pts,1)),triangles=t(tri))

All is left to do now is to plot the mesh to see what is our result. You can use the function shade3d from the package rgl:

shade3d(mesh,col="coral3",specular="orange")

 

Too Long, Didn’t read (or ‘just let me make some ammonites!’)

To make everything smooth and easy to use, what’s better than having all the above computation embedded into a handful of functions? The functions to produce the ammonite meshes are stored in the document DIYammonite_functions.R at this link. You can download them on your PC and then load them into your R environment using the function source and start having fun. The main function is diyAmmonite and it’s all you need to produce the mesh, given some values to the relevant parameters (a,b,turns and thick):

source("DIYammonite_functions.R")

#nspiral: the number of 3D points of which the spiral template is composed
#nsection: the number of 3D points for each circular section
#turns: the number of times the spiral changes direction while growing outward around the centre
#a: size parameter
#b: rate of radii increase (keep it between 0.1 and 0.3, or exceed this range and generate monsters)
#thick: reduces the radius of each section proportionally to its position along the spiral
#center: the x, y and z coordinates of the ammonite centre
ammo<-diyAmmonite(nspiral=200,nsection=30,turns=10,a=1.0,b=0.13,thick=2,center=c(0,0,0))
ammo3d<-ammo$mesh

shade3d(ammo3d,col="coral3",specular="orange")

 

Final considerations

Is this an accurate modelling of the ammonite shell morphology? - definitely not.

Can this produce an accurate picture of the ammonite morphological variability? - pretty sure not.

The functions here shown are better described as an experiment in coding and geometry, and the result of some lazy, rainy weekend hours. If you find any error/bug, please report it to me at .

I am sure this effort of mine is not in vain, and could be fruitfully used for educational purposes…or just to build your own collection of intact, almost life-like ammonite shells (possibly 3D-printable??), and admire them in all their beauty.