Counting tree rings

I saw this cross-section of a tree in La Connor, WA. Unusually, perhaps, for such a large tree, there was no nearby sign telling how old it was. Curious, I thought of a couple of ways to quickly get data so I could figure it out later.

Cross-section of a large tree

Image Analysis

The Google Camera app has a Panorama setting, which allowed me to get the above photo. So now it's possible for me to count the rings one by one, but that's boring, and isn't an easily repeatable solution.

Zoom of a cross-section of a large tree

If we look closer, though, each ring of the tree forms a convenient little ridge. At this lighting angle, they leave little shadows. We can work with this.

So I pulled out GIMP, cropped the image so that it included only the right-hand side of the tree, and fiddled with the Levels and Threshold settings to come up with a couple of black-and-white images.

Working on the tree in GIMP with levels

Working on the tree in GIMP with levels

Working on the tree in GIMP with threshold

Since the illumination and relative angle of the ridges changes across the diameter of the trunk, I also tried processing it in chunks.

Working on the tree in GIMP by taking regions

Working on the tree in GIMP with levels

The end result is a black and white image: a binary matrix.

Black and white view of the tree

Here's a zoom in on part of the above image.

Zoomed black and white view of the tree

I can now load this into R and find the number of tree rings like so:

library(png)
img<-readPNG("bw_tree_rings.png")

#Get only one of the dimensions. It doesn't matter 
#which one since this is a binary matrix so the RGB
#channels are all equal
img<-img[,,1]

#Check that a cross-section of the image is as expected
length(img[1,])

#If it wasn't, we could transpose. But it was, so we're good
#img<-t(img)

#Display the image
image(img, asp=1,axes=FALSE,useRaster=TRUE)

#Each row of the matrix is something like "0,0,0,1,1,1,0,1,1".
#Differentiating this gives "0,0,0,1,0,-1,0,1,0". Each value
#of "1" corresponds to a new ridge. Summing all the "1" values
#gives the number of tree rings. I do this for each row like so.
a<-apply(img, 1, function(x) sum(diff(x)>0))

#Print out a plot of tree rings versus image row
png('/z/tree-plot.png', width=6, height=6, units="in", res=300)
plot(a,xlab="Row", ylab="Tree Rings")
dev.off()

Which gives the following plot:

Tree rings by image row

Naturally, since not every row is a diamater, most of the above plot is rubish. But, looking at the original image, we see that the center of the tree is situated just off-center (vertically) of the image. Since the image is 2000 pixels tall, we are looking for a patch of data just to one side of the mid-point; here, that corresponds to the data points around row 950 or so.

This establishes the age of the tree as approximately 550–600 years old, though the latter is true only if you believe the highest data points.

Using the accelerometer

Since each tree ring formed a little ridge, another way to count the rings was to run AndroSensor at a high-frequency to collect my phone's accelerometer data. Holding the phone tightly in my hand to reduce damping, I ran my finger nail across a diameter of the tree. Each tree ring ridge resulted in a little bump, which the accelerometer would record.

The resulting data looks like this.

ACCELEROMETER X (m/s²) ACCELEROMETER Y (m/s²) ACCELEROMETER Z (m/s²) Time since start in ms YYYY-MO-DD HH-MI-SS_SSS
-1.152 6.942 6.2474 19 2015-08-16 18:24:28:327
-1.152 6.942 6.2474 36 2015-08-16 18:24:28:344
-1.152 6.942 6.2474 49 2015-08-16 18:24:28:357
-0.9712 6.9156 6.451 58 2015-08-16 18:24:28:366
-0.9712 6.9156 6.451 64 2015-08-16 18:24:28:372
-0.9712 6.9156 6.451 67 2015-08-16 18:24:28:375
-0.9712 6.9156 6.451 75 2015-08-16 18:24:28:383
-0.8442 7.1443 6.7132 78 2015-08-16 18:24:28:386
-0.8442 7.1443 6.7132 80 2015-08-16 18:24:28:388

Plotting the vector norm of the acceleration, the entire data set looks like this:

Accelerometer data for the tree

The first bit is me testing trying out the method. Then there's a quiet bit where I simply hold the phone with my finger resting on the tree. This gives me some calibration data so I can get a sense of how the accelerometer's background noise. Next, I begin running my fingernail across the tree trunk. There are occasional pauses as I relocate my body.

The calibration data looks a bit like this.

Accelerometer calibration data for the tree

The distribution of the calibration data about its mean is

Histogram of accelerometer calibration data for the tree

From the foregoing, it seems reasonable to classify every data point within 0.1 of the mean as being noise. Therefore, I manipulate the data by tossing out the negative side of the waveform (since every acceleration is followed by a deacceleration) and then removing all data beneath this noise threshold. I can then approximate the age by differentiating the signal and counting rises.

This gives 502 tree rings.

A hand count (I had to know!) gives 672 rings, so it seems that both of the methods I used were a bit losy.

The code for that is as follows:

library(ggplot2)

file='Sensor_record_20150816_182523_AndroSensor.csv'
dat<-read.csv(file, sep=';')
dat<-data.frame(dat)
colnames(dat)<-c('ax','ay','az','ms','date')
dat$at<-sqrt(dat$ax**2+dat$ay**2+dat$az**2)
head(dat)

#Show all the data
png('/z/tree-sound.png', width=6, height=3, units="in", res=300)
ggplot(dat)+geom_line(aes(x=ms,y=at))+ylab("Net Accel")+xlab("Time (ms)")
dev.off()

#Find the calibration period start and end times with
#plotting and some trial and error
calibration_period<-(7625<dat$ms & dat$ms<10295)
calib_dat<-dat[calibration_period,]
#plot(calib_dat$ms,calib_dat$at)

#Get the mean of the calibration data
datmean<-mean(calib_dat$at)

#Plot calibration period
png('/z/tree-calib.png', width=6, height=3, units="in", res=300)
ggplot(calib_dat)+
  geom_line(aes(x=ms,y=at))+
  xlab("Time (ms)")+
  ylab("Net Accel")+
  geom_hline(yintercept=datmean)
dev.off()

#Shift the data so it is centered about the mean
calib_dat$at<-calib_dat$at-datmean

#Plot a histogram of the shifted calibration data
png('/z/tree-hist.png', width=3,height=3,units="in",res=300)
hist(calib_dat$at, xlab="Deviation of net accel", ylab="Count", title="")
dev.off()

#Note that most of the noise falls within +/-0.1 of the mean, 
#and nearly all falls within +/-0.2

#Shift all of the data and take a peek
dat_period<-(10790<dat$ms & dat$ms<51489)
dat_dat<-dat[dat_period,]
dat_dat$at<-dat_dat$at-datmean

ggplot(dat_dat,aes(x=ms,y=at))+geom_line()+ylim(-1,1)#+xlim(25000,30000)

#Keep only the positive side of the waveform
dat_dat[dat_dat$at<0,]$at<-0

ggplot(dat_dat,aes(x=ms,y=at))+geom_line()+ylim(-1,1)+xlim(25000,30000)

#Remove the noise
dat_dat[dat_dat$at<0.1,]$at<-0

ggplot(dat_dat,aes(x=ms,y=at))+geom_line()+ylim(-1,1)+xlim(20000,30000)

#Count all the rises. This is an approximation of the age.
sum(diff(dat_dat$at)>0)

links