Tree Age from Image Analysis and Statistics

I found a nice snapshot of a cut tree that I took during a winter forest walk in my picture archive, and had the idea to automate the process of counting the tree rings in order to estimate the age of the tree when it was harvested.

Using VS Code, GitHub Copilot and some time to tinker and experiment, I wrote the some code to perform this analysis. The resulting Python program applies image processing techniques, detects circular features, and performs statistical analysis on the detected markers for such an estimation. It is by no means perfect – specifically finding the center of the rings proved challenging and error prone. The resulting estimates are, well, approximate.

The source code, input and process images, and additional samples to process can be downloaded here:

https://github.com/ferzkopp/TreeRingAnalysis

The program uses the following workflow:

Let’s walk through each step.

Since the source images may be small, the first step is an optional scaling operation to make sure it has  at least 3000 pixels in width. 

and to remove the noise some blurring is applied

Finding the center of the concentric rings – a task that is easy for a human looking at the image of a tree – is surprisingly difficult for an algorithm due to the nature of the image with noise, low contrast, vertical variations, and other inconsistencies. I tried initially several algorithms suggested by a quick Copilot query such as Hough Transform, Radial Symmetry Transform and Template Matching before settling on a gradient decent search using a cost function that is designed to find the center of concentric structures (like tree rings) based on edge data in the image.

# Cost function to evaluate the center

def cost_function(center, edges):

    x_center, y_center = center

    y_indices, x_indices = np.where(edges > 0)  # Get edge points

    radii = np.sqrt((x_indices - x_center)**2 + (y_indices - y_center)**2)

    mean_radius = np.mean(radii)

    cost = np.sum((radii - mean_radius)**2)  # Minimize the deviation from the mean radius

    return cost

# Gradient descent to find the best center

def find_best_center(edges, initial_center):

    result = minimize(

        cost_function,

        initial_center,

        args=(edges,),

        method='Powell'

    )

    return result.x

The input to this step is a “binarized” image that shows the rings around the center as bright pixels. The following processing was performed to create such a black and white image:

• cut out the center quarter of the image, since it is most likely to contain the center of the rings
  
• convert the color image into HSV (Hue, Saturation, Value) format
  calculate a histogram of the V component (brightness)
• binarize values using a percentile threshold (values below X^th-percentile become white pixels)
• mask out all parts not within a centered circle (black pixel)

Here is an example of the binarized version of the above image source:

and this is the result of the center search based on the masked version of this image:

The green dot indicates the initial estimate at the start of the search (image center) and the red dot indicates the final location (estimate of tree ring center). As one can see, this estimate is far from perfect, but sufficiently accurate to proceed.

Now that we have a center, the next step of the algorithm can be applied: creating time-series of brightness values along rays emanating from the center. First a grayscale version of the source image is created and enhanced using histogram equalization.

Then the image is sampled along 360 rays projecting from the center in 1 degree intervals to create 360 individual brightness series.

Here is a close-up showing the dark (orange dots) and bright (blue dots) areas detected in the image.

Each brightness series is then processed as follows:

• the series is frequency filtered with a low-pass filter to remove high frequency variations
• the filter is using a dynamic cutoff, filtering more close to the center (which has broader spaced rings) and less in the outer parts (which usually has tighter ring spacings) using a linear cutoff response
  
• then on to detect peaks and troughs in the resulting brightness curve, treating each peak (or trough) as a possible ring,
• and finally by collecting the count of the peaks (or troughs) in each series as an estimate along that particular ray

fft_values = np.fft.fft(brightness_values)
frequencies = np.fft.fftfreq(len(brightness_values))
dynamic_cutoff = cutoff + (cutoff * 2) * (frequencies / max(frequencies))
fft_values[np.abs(frequencies) > dynamic_cutoff] = 0
smoothed_brightness_values = np.fft.ifft(fft_values).real
troughs, _ = find_peaks(-smoothed_brightness_values, prominence=prominence)
peaks, _ = find_peaks(smoothed_brightness_values, prominence=prominence)

The following chart shows the resulting peak and trough detection results for 3 rays (0/45/90 degrees):

The final step is to statistically analyze the set of peak and trough counts (collected as markers) that were found, 720 in total, to find a possible range for our tree age using a mean, mode and maxima in the set.

This statistical averaging “smooths” the imperfections in the algorithms.

The program outputs the final estimate:

...

The mean of the marker series is 55.
The mode of the marker series is 58 (with a count of 75).
The maximum value of the marker series is 67.
The estimated age of the tree is between 55 and 67 years.

I counted the tree rings in the image by hand and estimate the actual age of the tree to be about 72 years, so the algorithm was underestimating. However, even manually determining the ring count from the image is challenging and my own estimate may not be entirely correct.
Close enough for me. 😉