Isolating phyllotactic patterns embedded in the secondary growth of sweet cherry (Prunus avium L.) using magnetic resonance imaging

Isolating phyllotactic patterns embedded in the secondary growth of sweet cherry (Prunus avium L.) using magnetic resonance imaging

[    [    [    [    [ \orgnameDepartment of Computational Mathematics, Science and Engineering, Michigan State University, \cityEast Lansing, MI, \cnyUSA \postcode48823 \orgnameDepartment of Computational Mathematics, Science and Engineering, Michigan State University, \cityEast Lansing, MI, \cnyUSA \postcode48823 \orgnameDepartment of Mathematics, Michigan State University, \cityEast Lansing, MI, \cnyUSA \postcode48823

Background: Epicormic branches arise from dormant buds patterned during the growth of previous years. Dormant epicormic buds remain on the surface of trees, pushed outward from the pith during secondary growth, but maintaining vascular connections. Epicormic buds can be reactivated, either through natural processes or intentionally, to rejuvenate orchards and control tree architecture. Because epicormic structures are embedded within secondary growth, tomographic approaches are a useful method to study them and understand their development.

Results: We apply techniques from image processing to determine the locations of epicormic vascular traces embedded within secondary growth of sweet cherry (Prunus avium L.), revealing the juvenile phyllotactic pattern in the trunk of an adult tree. Techniques include breadth-first search to find the pith of the tree, edge detection to approximate the radius, and a conversion to polar coordinates to threshold and segment phyllotactic features. Intensity values from Magnetic Resonance Imaging (MRI) of the trunk are projected onto the surface of a perfect cylinder to find the locations of traces in the “boundary image”. Mathematical phyllotaxy provides a means to capture the patterns in the boundary image by modeling phyllotactic parameters. Our cherry tree specimen has the conspicuous parastichy pair , phyllotactic fraction 2/5, and divergence angle of approximately degrees.

Conclusions: The methods described not only provide a framework to study phyllotaxy, but for image processing of volumetric image data in plants. Our results have practical implications for orchard rejuvenation and directed approaches to influence tree architecture. The study of epicormic structures, which are hidden within secondary growth, using tomographic methods also opens the possibility of studying the genetic and environmental basis of such structures.



addressref=cmse, ]\initsME\fnmMitchell \snmEithun addressref=hort, ]\initsJL\fnmJames \snmLarson addressref=hort, ]\initsGL\fnmGregory \snmLang addressref=hort,cmse, ]\initsDHC\fnmDaniel H. \snmChitwood addressref=cmse,math, ]\initsEM\fnmElizabeth \snmMunch


Magnetic resonance imaging \kwdSweet cherry \kwdPhyllotaxy \kwdParastichy \kwdThresholding

1 Introduction

Plants increase in length from apical meristems during primary growth. Located at the shoot tip, the shoot apical meristem is the site of cell division that produces leaf and bud primordia. Newly divided cells elongate, pushing the apical meristem upward. As the shoot continues to grow, leaf primordia cells divide, differentiate, and expand into leaves that subtend axillary buds. The point at which each leaf is attached to the shoot constitutes a node. In sweet cherry (Prunus avium L.), leaves form at each node in the year that a shoot develops. In non-juvenile plants, the axillary buds at nodes at the base of the shoot may begin differentiating into solitary flower buds. These bloom, develop fruit if pollinated and fertilized, and upon abscission of the fruit, the node becomes void of apparent vegetative or reproductive buds for subsequent growth. These are called “blind” nodes. The remaining non-basal, majority of nodes on the new shoot form a single vegetative bud at each node. In spring of the year after the shoot’s formation, each of these buds usually produces five to eight leaves that arise from very closely-packed nodes, producing a spur (a short, modified branch) along the rest of the length of the shoot, except for the terminal shoot apical meristem that again elongates to form a new section of shoot. Buds on some of the uppermost nodes of the original shoot also may elongate into new lateral shoots, rather than form spurs [1]. In the orchard, trees are considered mature once they have filled their allotted orchard space and have all of their reproductive components; in modern, high density plantings, this is three to five years.

The shoot apical meristem forms nodes in spiral patterns resembling cylindrical helices. Phyllotaxy is the study of the arrangement of plant organs during their development. These organs can be branches, leaves, vascular traces, or other features associated with nodes. Frequently, organ primordia form spiral patterns called parastichies, which can be characterized by Fibonacci-type sequences of numbers. For centuries, mathematicians have been interested in building models to describe the geometry of phyllotactic patterns [2]. In these models, phyllotactic parameters are numbers that capture the layout and spacing of primordia characteristic of different plant species and their stages of development [3].

Phyllotactic patterns are often studied in two forms: centric and cylindrical [2]. In the centric representation, parastichy spirals emanate from a central point, like the capitulum of a sunflower. In the cylindrical model, parastichies are helices on the surface of a cylinder, like in pineapples or pinecones. In the cylindrical representation, it is convenient to unwrap the surface of the cylinder and view the primordia as a cylindrical lattice in the plane, called a Bravais lattice, in which intersections of parastichies are primordia.

After a shoot has elongated and the phyllotaxy of the nodes patterned, woody plants undergo secondary growth to increase in girth. Axillary buds usually remain dormant in the year of initiation, but when the primary shoot is damaged or is growing extremely rapidly, the axillary bud may elongate into a new lateral shoot. Axillary buds that remain dormant will eventually become engulfed by secondary growth of the stem and persist beneath the bark as an epicormic bud meristem [4, 5]. Epicormic bud meristem cells divide and elongate with radial growth, maintaining their presence just beneath the bark, and leaving a vascular trace back to the pith, shown in Fig. 2. Epicormic traces are 2-5 mm high [6] and occur perpendicular to the pith [4, 7]. Epicormic buds may sprout into epicormic branches following a stress such as fire [8], insect defoliation [9], wind damage [10], competition [11], or pruning [12]. If the primary epicormic meristem becomes damaged, it may split via the initiation of subtending secondary bud meristems. Epicormic buds can be reactivated under the right conditions, and serve as a reserve of future potential shoot growth that can be used to rejuvenate orchards by farmers.

Image processing techniques have previously been used for dendochronology, the study of tree rings and their features [13, 14] and to manually locate rameal traces in oak trunks [7]. Semiautomatic feature identification on trees using image processing has been used for tree ring identification using the Sobel operator [15]. In this work, we isolate the patterning of epicormic traces embedded in secondary growth from a magnetic resonance imaging (MRI) scan of an eight year old sweet cherry tree. An image processing framework is implemented that finds the pith and radius of the trunk across slices. Using these dimensions, a polar coordinate conversion of each slice reveals x-ray dense regions corresponding to epicormic traces. A blob detection algorithm segments the epicormic features and the resulting phyllotactic parameters are estimated. Our work reveals the juvenile phyllotactic pattern of a sweet cherry tree embedded within 8 years of secondary growth. The study of epicormic features has implications for orchard rejuvenation, and the analysis methods presented provide an empirically based method to measure phyllotactic patterning and isolate anatomical features in plants.

2 Materials and Methods

2.1 Plant material and Magnetic Resonance Imaging (MRI)

At bloom (late April) in 2016, the top of an 8-year-old sweet cherry tree, cultivar ’NY 119’, planted at Michigan State University’s Clarksville Research Center and trained to a standard central leader canopy architecture, was removed with a chain saw 1.5 m from the ground. All lateral branches below that point were removed to promote the sprouting of epicormic buds along the remaining trunk length. Prior to and following this topping, the tree was managed with standard fertility, irrigation, and pest management procedures with the rest of the ha experimental orchard. In August, a 1.14 m-long section of trunk was fully removed just above the graft union and dried at room temperature in the laboratory. In December, this trunk section was scanned at slice thickness set of 0.625 mm, at Michigan State University’s Department of Radiology (East Lansing, MI) using a whole body magnetic resonance scanner (GE Signa HDX 3.0T, Chicago, IL). The scan took slightly over a minute and produced 1871 images, each of which is pixels. Contrast in the scanned image is created by differences in moisture content; contrast is greater in air-dried than fresh wood [16]. The specimen and the scan are shown in Fig. 1. Table 1 gives estimates for the size and shape of the specimen.

Measurement Value
# MRI slices 1871
tallest height 115.13 cm
shortest height 110.31 cm
median radius 4.72 cm 0.01 cm
median circumference 29.64 cm 0.08 cm
Table 1: Sweet cherry specimen statistics. Uncertainty is the standard error of the median.

2.2 Image Processing

MRI generates a series of dicom files corresponding to slices of the object being scanned. Each slice was taken parallel to the the ground so that each is an image of rings of the tree at a fixed height. For ease of notation, we rescale the intensity values over all pixel arrays in a scan so that they span .

Mathematically, a pixel slice (or “slice”) is a function

This function maps pixel positions to intensity value. For example, if the the pixel in the upper-left hand corner of the slice has intensity value , then we write . Notice that the notion of a pixel slice is equivalent to a 2D image, making it amenable to methods from the field of image processing.

2.2.1 Thresholding

Thresholding is a fundamental image processing technique which generates a binary image intensity threshold. Given a pixel slice , select a threshold . Then define a thresholded slice by

Thresholding can be used to segment an image or reveal features of interest by eliminating pixels whose value is below . There are many ways to select a threshold. We use Otsu’s method, which assumes an underlying distribution of intensity values from two classes and iterates to minimize intra-class variance [17].

2.2.2 Edge Detection

Given an input image, edge detection algorithms generate a new image marking probable edges. The Sobel edge detection algorithm works by convolving an image with the Sobel operator to approximate the image gradient, which highlights areas of abrupt change in intensity [18]. Previously the Sobel operator has also been used for tree ring identification in images of cross-sections [15]. We apply the Sobel algorithm to estimate the radii of slices in Section 3.3.

2.2.3 Graphs and Breadth-First Search

In order to traverse pixel locations in a pixel slice, we also think of each slice as a mathematical graph. A graph is a mathematical construct consisting of a set of objects and connections between them. The objects are called edges and the connections are called vertices. A graph is usually denoted by , where is the set of vertices and is the set of edges.

To represent a slice as a graph, let be all coordinate pair inputs (pixels). An edge exists between two coordinate pairs if exactly one entry in the pairs differs by exactly one. In other words, there exists an edge from a pixel to its neighbors in the up, down, left and right directions. More formally, let be an pixel slice and define a set of directions . Then define the slice graph by

The restrictions and on elements of the vertex set ensure that there actually exists a pixel neighboring in the direction given by .

Using a slice graph to represent a pixel slice allows us to traverse the image and look at pixels in a certain order. In particular, in part of our procedure we use breadth-first search (BFS), which we will summarize here. For a more detailed explanation of BFS in image processing, see [19].

Search algorithms for graphs differ by the order in which we look at pixels. BFS is a search algorithm for graphs which looks at nearby neighbors before considering neighbors that are farther away (“breadth-first” is in opposition to “depth-first”). In practice, on a slice graph, this means that if we are seeking a pixel that meets some condition nearby a seed , we first look at the the intensity of its neighbors. After considering each neighbor of we look at the neighbors of the neighbors of . Depending on the application, this process continues until we find a pixel the meets the condition or we find the largest region containing that meets the condition. This technique will be used in Sec. 2.2 to determine the centroid of the pith in each slice.

2.3 Mathematical Phyllotaxy

The goal of mathematical phyllotaxy is to describe emergent spiral and other patterns of lateral organs. For the cylindrical representation, this includes viewing phyllotactic patterns as helices on a cylinder or sets of parallel lines. The idealized geometric model, describing the placement of the primordia and parastichies is described by the Fundamental Theorem of Phyllotaxy [20, 21], from which Jean’s pattern determination table can be used to estimate phyllotactic parameters; e.g., using an allometry-based model [22]. Fibonacci numbers associated with parastichies can be seen as fixed points in dynamical systems derived from optimal packing assumptions [23]. For a full review of links between mathematical and molecular phyllotaxy, see [24].

Besides their ubiquity in developing plant structures, phyllotactic features appear in other natural phenomena such as capillary structures during evaporation [25] and the cylindrical phyllotaxy of carbon nanontubes [26]. For a discussion of the universality of Fibonacci patterns in nature, see [27]. In what follows we define the several phyllotatic parameters used to summarize the arrangement of plant organs.

2.3.1 Parastichy

Adopting a convention from [2], we call the plant organs that comprise phyllotactic patterns primordia. In the case of our work, the primordia are epicormic buds near the epidermis connected by radial, vascular traces to the pith.

A parastichy is a set of primordia that form a spiral arm. All parastichies that run in the same direction form a family and a parastichy in a family of spirals is an -parastichy. Fig. 3 shows families of 5- and 8-parastichies in synthetic data. The 1-parastichy is called the genetic spiral and contains all primordia. A contact parastichy is one that is derived from the shape of the primordia. For example, each hexagonal primordia on the outside surface of a pineapple suggests three contact parastichies.

A parastichy pair is formed by parastichies in one direction and parastichies in the other. Parastichy pairs are often consecutive numbers in a Fibonacci-type sequence. Define a sequence by , and , where and . In normal phyllotaxy parastichy numbers come from the sequence and Fibonacci numbers arise in the case that and [2]. We restrict ourselves to the Fibonacci sequence .

A visibly opposed parastichy pair intersects only at primordia. There are many of these pairs in a phenomenon known as “rising phyllotaxy”. In order to choose a definite parastichy pair, we are interested in the conspicuous parastichy pair, which is a visibly opposed parastichy pair such that the angle of intersection between the two families is closest to [2]. This pair is not necessarily unique, but it is the most “conscipuous” in the sense that the corresponding families of parastichies are the most noticeable.

2.3.2 Parameters

Since the cherry tree trunk is cylindrical, we limit our discussion of phyllotactic parameters to cylindrical phyllotaxy. In the cylindrical representation the coordinates of the primordia on the surface of the cylinder are given by

where for (this notation implies an arbitrary position for the ray ). The divergence angle is defined as the average difference between successive angle measurements, i.e.

For clarity we also define the divergence fraction as , which is the fraction of the angular breadth of the arc made by the divergence angle (in some literature this is the definition of the divergence angle [2]).

Another useful parameter is the vertical distance between successive primordia, or the rise, defined by the internode distance, . The conspicuous parastichy pairs are a function of both the divergence angle and the rise. Hence, both parameters are necessary to fully to describe the phyllotaxy of the system.

Using the divergence angle , spiral nodes can be generated by and cylindrical lattice points can be generated by . Examples of both centric and cylindrical phyllotaxy are shown in Fig. 3.

2.3.3 Phyllotactic Fractions

In contrast to paristichy pairs, a more traditional way to describe a cylindrical phyllotactic system is to use a phyllotactic fraction to approximate the angular differences between nodes [28]. A phyllotactic fraction (sometimes called “the phyllotaxis” of a species) is determined by the number of turns of a spiral, , of successive leaves to reach the th node directly above the starting node. In the fraction, and are always every other number in a Fibonacci-type sequence from normal phyllotaxy. For example, cherry’s phyllotactic fraction, derived from the Fibonacci sequence, is ; this means that a node will have a node form above it after approximately two spirals of five successive nodes. A cherry stem will then have five vertical ranks of nodes called orthostichies [29]. The phyllotactic fraction is meant to approximate the divergence angle. For example, the phyllotactic fraction 2/5 implies a divergence angle of .

2.3.4 Fundamental Theorem of Phyllotaxy

Using the divergence fraction, we can identify parastichies in an idealized model describing the pattern of primordia.

First, creftype 1 helps identify which primordia belong to which parastichies. After sorting the primordia by their radial angle and then labelling the primordia by their position, we call two primordia adjacent if they are next to each other in this list.

Theorem 1 (The Bravais-Bravais Theorem, 1837 [2]).

On an -parastichy of a phyllotactic spiral pattern, the numbers of two adjacent primordia differ by .

It follows from creftype 1, for example, that the parastichies in the 2-parastichy family will have labels and . To identify which parastichies are visibly opposed is characterized by creftype 2. For ease of notation let denote the nearest integer to the real number .

Theorem 2 (Fundamental Theorem of Phyllotaxy, 1994 [2]; Revised 2012 [21]).

The following are equivalent:

  1. The parastichy pair is visible and opposed, where and .

    1. , , and , or

    2. , , , , or

    3. , , , , or

    4. , , , where are the unique integers , and such that .

Notice that bounding the divergence angle in a particular way is necessary and sufficient to determine if a parastichy pair is visible and opposed. The following theorem restates the creftype 2 in terms of the Fibonacci sequence.

Theorem 3 (Adler’s Theorem, 1974 [30]).

Let denote the Fibonacci sequence and define an interval

The parastichy pair is visible and opposed if and only if .

Note that , where is the golden ratio. Hence, if is visble and opposed for all , then . The golden angle is defined as .

3 Results

Using image processing techniques and mathematical phyllotaxy, we propose a method to algorithmically determine the conspicuous parastichy pair found in cylindrical phyllotaxy. This method could be the basis of an automatic method to retrieve the locations of primordia in 3D scans of plants with cylindrical phyllotaxy, given that the primordia are revealed by the distribution of intensity values. Again let denote the nearest integer to the real number .

3.1 Pre-Processing

The cherry tree scan is comprised of 1871 dicom files, each of which contains a pixel slice. Each slice represents a thickness of 0.625 mm and each pixel a 0.351562 mm 0.351562 mm square. Hence, each voxel cube in the image has a volume of about 0.077 mm. Using the pydicom Python package [31], we converted dicom files to numpy matrices [32].

3.2 Pith Finding Algorithm

Since the tree trunk is not a perfect cylinder, the location of the tree’s pith in each slice is different. To be able to orient the object in a coordinate system, we first identify the location of the pith centroid in each slice which we will use later to convert each slice to polar coordinates.

To keep the description of procedure as general as possible, suppose we have a scan with slices denoted by . Given the center of the pith in , we iteratively locate the center of the pith in slices using a search algorithm. First define a range of intensity values which represents the intensity values found in the pith in the whole scan, since the pith is a region (collection of coordinate pairs) and just one coordinate pairs. In our image, we chose . Assuming there exists a boundary around the target pith region in each slice whose intensity values are not in , we can use to search for the pith region.

We require user input to determine a starting location for the pith in the first slice, . As this will likely not be exactly the centroid of that pith, we will call this initialized input . Let be the center of the pith in slice . Using BFS in slice , we spread outward from until we find a coordinate pair such that . Note that it is possible for to be the starting coordinate if . This puts us inside the pith region of slice . We find the whole pith region by using another BFS to find the entire connected region containing , such that all each coordinate pairs in the region have intensity values in . Finally, the centroid of this region is marked as the center of the pith for slice . The full algorithm is outlined in Algorithm 1.

Data: Slices , initial center , intensity range
Result: Center coordinates
for i=1,…,N do
       Using BFS, choose
       Using BFS, find , the largest connected region of such that and for all
       Set to be the number of pixels in
       Set and
end for
Algorithm 1 Algorithm for finding pith centroid locations

3.3 Radius Estimation

To estimate the average radius of the cherry tree we identify the edge of the tree in each slice and then compute the distance to the pith location. This process is outlined in Fig. 4. Specifically we apply the Sobel edge detection algorithm to identify probable edges. Since the edge detection algorithm also identifies rings and traces in the interior of the tree, we threshold the edges and keep only values above the 99th percentile, leaving values on the boundary of the tree. Then we compute the median distance from these points to the pith centroid to estimate the radius of a slice. Let denote the estimated radius of slice . We estimate the overall radius of the tree as the median of over all .

3.4 Polar Conversion

Using the coordinate pairs for pith locations produced by Algorithm 1 and the estimated radii, we convert each slice to polar coordinates. Then the radial, wedged-shaped traces in the slices become vertical blocks and are easier to identify.

Specfically, for each slice , define a polar slice by

for , where are small positive numbers. Selected pixel slices and polar slices are shown in Fig. 5.

3.5 Boundary Image

Next we construct a radial summary of the polar slices that we call the boundary image. Notice that the traces in Fig. 5. have higher intensity values than the other parts of the tree. This means the intensity values in a column will tend to be higher if a trace is present. Hence, a summary statistic of the intensity values in each column may capture whether or not a trace exists there.

For each polar image , define an intensity range that denotes the foreground of the image. Define the boundary image by

Hence, the mean pixel intensity along columns in polar slices become rows in the boundary image (the intensity range restricts us to consider only pixels that are a part of the tree and not pixels that are in the background). Note that other summary statistics such as the median or the percentile may be used in place of the mean.

The boundary image is shown in Fig. 6(A). The boundary rays and glue together so that the boundary image is the boundary of a cylinder representing the cherry tree. Notice the phyllotactic patterns that emerge in the lattice structure of the image.

3.6 Blob Detection

Each blob in the boundary image corresponds to an epicormic trace in the cherry tree. To identify the centroid of each blob, we threshold the image using Otsu’s method. Then we compute the centroid of each region to find the trace locations, shown in Fig. 6(B). Notice that one blob (colored in red) was too small to be recognized. The centroid of this blob was found by using a less restrictive threshold and was added to the final dataset.

3.7 Phyllotaxy Parameter Determination

Using the centroids of the blobs in Fig. 6(B), we estimate phyllotactic parameters in the context of both the Fundamental Theorem of Phyllotaxy (creftype 2) and phyllotactic fractions. Start by sorting and labeling the 36 nodes by their location so that node 0 is near the bottom of the trunk and node 35 is near the top. The divergence angle, estimated by the mean change in angle between successive primordia is , which gives a divergence fraction of 0.397 0.003. To see the distribution of local divergence angles and rise, see Fig. 7. Assume that is the actual divergence fraction of the system. Since , it follows from Adler’s Theorem (creftype 3) that are all visible and opposed parastichy pairs. Also, 1, 2, 3, 5 and 8 are parastichy numbers in the sense that are part of visibly opposed parastichy pairs for the specimen. However, the standard error on suggests that the true mean could be greater than , which would exclude the 8-parastichies from consideration. Fig. 8 shows all parastichies in 3D and on the Bravais lattice.

Since the second number in a phyllotactic fractions corresponds to a parastichy number, 1/2, 1/3, 2/5 and 3/8 are all phyllotactic fractions of the system (see [2, §2.2.2]). Traditionally, the phyllotactic fraction , chosen to represent the phyllotaxy of a cylindrical system, is derived from the following observation: take two nodes whose labels differ by and note that they are approximately above each other (having the same values) by going exactly times around the genetic spiral. In this sense, the “phyllotaxy of the system” for the cherry specimen can be described as 2/5. We propose a more precise way of determining a single phyllotactic fraction: take the Fibonacci fraction closest to the divergence fraction . In other words, for phyllotaxy based on the Fibonacci sequence, let be the conspicuous phyllotactic fraction, where . In our case, , and since 2/5 is closest to .

To determine the conspicuous parastichy pair, we calculate angles between parastichies. Using the Bravais-Bravais Theorem (creftype 1), draw -parastichies and -parastichies using linear interpolation and measure the angles at each intersection. The pair of and for which this angle is closest to is the conspicuous parastichy pair (all parastichy families are shown in Fig. 8). Using the law of cosines, we compute the average intersection angles between families of parastichy pairs and convert each to the “small angle”: for (2,3); for (3,5); and for (5,8). Thus, by definition, the parastichy conspicuous parastichy pair is . The 2-family and 3-family of parastichies are shown in the Bravais lattice in Fig. 6(C). A summary of phyllotactic parameters found for this sample are shown in Table 2.

Parameter Value
number of primordia 36
handedness counterclockwise
divergence angle () (2.495 rad 0.021 rad)
divergence fraction () 0.397 0.003
average rise () 2.896 cm 0.182 cm
visible & opposed parastichy pairs (1,2), (2,3), (3,5), (5,8)
parastichy numbers 1, 2, 3, 5, 8
conspicuous parastichy pair (2,3)
conspicuous phyllotactic fraction
Table 2: Phyllotactic parameters with standard error. The “handedness” of the system is the direction of the genetic spiral going from the bottom of the sample. The families of parastichies alternate between clockwise and counterclockwise direction.

4 Discussion

The 36 nodes detected in our sample were patterned in a single year. The traces connecting to the epicormic buds traverse eight years of secondary growth. In the orchard, trees are planted as a “whip” of a genetically compound tree, that is, a one-year-old nursery tree is comprised of the shoot of the scion genotype (fruiting variety) that grew in the nursery the previous year from a bud that was grafted onto a rootstock genotype. A whip nursery tree is typically 1 to 1.5 m tall. In this experiment, the number (8) of annual growth rings at the top of the scanned trunk section matched the number at the bottom; therefore, each node in the section of the scanned trunk originated in the same year, all on the original whip nursery tree. To our knowledge, this is the first accounting of a single growing season’s complement of sweet cherry nodes from origin through eight years of trunk growth that also identified the persistence of epicormic bud traces. That is, the phyllotactic patterning that occurred during juvenile shoot growth remains in a mature tree, and can be detected using a combination of magnetic resonance imaging and image processing approaches.

4.1 Conclusions

Our results have practical implications for orchard rejuvenation and directed approaches to influence tree architecture. As orchards age, yield and fruit quality can begin to decline if fruiting sites become shaded and/or portions of the tree become infected by diseases. Epicormic buds serve as a “bank” for new branches to sprout and rejuvenate orchards. Magnetic resonance imaging provides a clear picture of how full that bank was after eight years. This could facilitate study of how persistence of epicormic buds may be affected by cultivar, vigor, harsh winters, drought, disease, and spring freezes. Determination of phyllotaxis provides growers with the potential to identify where an epicormic bud is located so that orchard training measures may be attempted to force a branch to sprout at that location. This could help growers fill “holes” in tree canopies to increase fruiting potential and lengthen the life of the orchard. The study of epicormic structures, which are hidden within secondary growth, using tomographic methods also opens the possibility of studying the genetic and environmental basis of such structures.

Our results also provide an empirical way to measure phyllotactic parameters. Much work in phyllotaxy has focused on generative models. But the image processing techniques presented here provide a method to isolate nodes and place them in the context of shoot growth using anatomical features. This is even possible if the features are difficult to discern by eye, or embedded internally within extensive secondary growth. With larger sets of node locations, Fourier Methods could further automate the process of finding parastichy numbers [33, 34]. Isolating features in mature specimens can lead to insights regarding the developmental history of a plant, which is crucial to manipulate plant forms in a directed way and mechanistically understand the origins of morphology. Automated methodology to model imaging features—–from MRI or otherwise—–into a developmental context is an important first step towards an empirical mathematical framework for measuring plant morphology.


Author’s contributions

ME, JL, GL, DHC, and EM designed the project and contributed to writing and editing the manuscript. JL and GL collected plant materials and performed the MRI scan. ME, DHC, and EM analyzed data. ME and EM developed new analysis tools and performed mathematical analyses.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

Codes are available on Github (, and raw data are available on the figshare repository (


This project was supported by the USDA National Institute of Food and Agriculture, and by Michigan State University AgBioResearch. The work of EM was supported in part by NSF grants DMS-1800446 and CMMI-1800466.


  • [1] Quero-Garcia, J., Iezzoni, A., Pulawska, J., Lang, G.A.: Cherries: Botany, Production and Uses. CABI. Google-Books-ID: 5xsxDwAAQBAJ (2017)
  • [2] Jean, R.V.: Phyllotaxis by Roger V. Jean (1994). doi:10.1017/CBO9780511666933. /core/books/phyllotaxis/272D9010BE175D26B61D5A2ED8D87A3C Accessed 2018-07-02
  • [3] I. Vakarelov, I.: Method for practical assessment of cylindrically represented spiral phyllotaxis (2018)
  • [4] Busgen, M., Munch, E.: The Structure and Life of Forest Trees. London, Chapman & Hall (1929). OCLC: 704188260
  • [5] Stone, E.L., Stone, M.H.: ”Dormant” Versus ”Adventitious” Buds. Science 98(2533), 62–62 (1943). doi:10.1126/science.98.2533.62. Accessed 2018-10-31
  • [6] Fontaine, F., Druelle, J.-L., Clément, C., Burrus, M., Audran, J.-C.: Ontogeny of proventitious epicormic buds in Quercus petraea. I. In the 5 years following initiation. Trees 13(1), 54–62 (1998). doi:10.1007/PL00009737. Accessed 2018-10-31
  • [7] Colin, F., Mothe, F., Freyburger, C., Morisset, J.-B., Leban, J.-M., Fontaine, F.: Tracking rameal traces in sessile oak trunks with X-ray computer tomography: biological bases, preliminary results and perspectives. Trees 24, 953–967 (2010). doi:10.1007/s00468-010-0466-1
  • [8] Burrows, G.: Syncarpia and Tristaniopsis (Myrtaceae) possess specialised fire-resistant epicormic structures. Australian Journal of Botany 56(3), 254–264 (2008). doi:10.1071/BT07164. Accessed 2018-10-31
  • [9] Piene, H., Eveleigh, E.S.: Spruce budworm defoliation in young balsam fir: the “green” tree phenomenon (1996). Accessed 2018-10-31
  • [10] Cooper-Ellis, S., Foster, D.R., Carlton, G., Lezberg, A.: Forest Response to Catastrophic Wind: Results from an Experimental Hurricane. Ecology 80(8), 2683–2696 (1999). doi:10.2307/177250. Accessed 2018-10-31
  • [11] Nicolini, E., Chanson, B., Bonne, F.: Stem Growth and Epicormic Branch Formation in Understorey Beech Trees ( Fagus sylvatica L.). Annals of Botany 87(6), 737–750 (2001). doi:10.1006/anbo.2001.1398. Accessed 2018-10-31
  • [12] O’Hara, K.L., Berrill, J.-P.: Epicormic sprout development in pruned coast redwood: pruning severity, genotype, and sprouting characteristics. Annals of Forest Science 66(4), 409–409 (2009). doi:10.1051/forest/2009015. Accessed 2018-10-31
  • [13] Levanic, T.: Atrics – A New System for Image Acquisition in Dendrochronology. Tree-Ring Research 63, 117–122 (2009). doi:10.3959/1536-1098-63.2.117
  • [14] Sundari, P.M., Kumar, S.B.R., Sahayaraj, A.J.: An Approach for Analyzing the Factors Recorded in the Tree Rings Using Image Processing Techniques. In: 2017 World Congress on Computing and Communication Technologies (WCCCT), pp. 236–239 (2017). doi:10.1109/WCCCT.2016.65
  • [15] Henke, M., Sloboda, B.: Semiautomatic tree ring segmentation using Active Contours and an optimised gradient operator. Forestry Journal 60 (2014). doi:10.2478/forj-2014-0020
  • [16] Freyburger, C., Longuetaud, F., Mothe, F., Constant, T., Leban, J.-M.: Measuring wood density by means of X-ray computer tomography. Annals of Forest Science 66(8), 804–804 (2009). doi:10.1051/forest/2009071. Accessed 2018-10-31
  • [17] Otsu, N.: A Threshold Selection Method from Gray-Level Histograms. IEEE Transactions on Systems, Man, and Cybernetics 9(1), 62–66 (1979). doi:10.1109/TSMC.1979.4310076
  • [18] Gupta, S., Mazumdar, S.G.: Sobel Edge Detection Algorithm 2(2), 6 (2013)
  • [19] Silvela, J., Portillo, J.: Breadth-first search and its application to image processing problems. IEEE Transactions on Image Processing 10(8), 1194–1199 (2001). doi:10.1109/83.935035
  • [20] Jean, R.V.: A mathematical model and a method for the practical assessment of the phyllotactic patterns. Journal of Theoretical Biology 129(1), 69–90 (1987). doi:10.1016/S0022-5193(87)80204-7. Accessed 2018-07-02
  • [21] Swinton, J.: The Fundamental Theorem of Phyllotaxis revisited. arXiv:1201.1641 [q-bio] (2012). arXiv: 1201.1641. Accessed 2018-10-31
  • [22] Jean, R.V.: Allometric relations in plant growth. Journal of Mathematical Biology 18(3), 189–200 (1983). doi:10.1007/BF00276086. Accessed 2018-07-31
  • [23] Atela, Golé, Hotton: A Dynamical System for Plant Pattern Formation: A Rigorous Analysis. Journal of Nonlinear Science 12(6), 641–676 (2003). doi:10.1007/s00332-002-0513-1. Accessed 2018-09-30
  • [24] Pennybacker, M.F., Shipman, P.D., Newell, A.C.: Phyllotaxis: Some progress, but a story far from over. Physica D: Nonlinear Phenomena 306, 48–81 (2015). doi:10.1016/j.physd.2015.05.003. Accessed 2018-09-30
  • [25] Chen, C., Duru, P., Joseph, P., Geoffroy, S., Prat, M.: Control of evaporation by geometry in capillary structures. From confined pillar arrays in a gap radial gradient to phyllotaxy-inspired geometry. Scientific Reports 7(1), 15110 (2017). doi:10.1038/s41598-017-14529-z. Accessed 2018-09-29
  • [26] Fan, X.-D., Bursill, L.A.: Principles for structure analysis of carbon nanotubes by high-resolution transmission electron microscopy. Philosophical Magazine A 72(1), 139–159 (1995). doi:10.1080/01418619508239587. Accessed 2018-06-28
  • [27] Shipman, P.D., Sun, Z., Pennybacker, M., Newell, A.C.: How universal are Fibonacci patterns? The European Physical Journal D 62(1), 5–17 (2011). doi:10.1140/epjd/e2010-00271-8. Accessed 2018-10-01
  • [28] Allard, H.A.: Some Aspects of the Phyllotaxy of Tobacco. Journal of Agricultural Research 64(1), 49–55 (1942)
  • [29] Okabe, T.: Biophysical optimality of the golden angle in phyllotaxis. Scientific Reports 5, 15358 (2015)
  • [30] Adler, I.: A model of contact pressure in phyllotaxis. Journal of Theoretical Biology 45(1), 1–79 (1974). doi:10.1016/0022-5193(74)90043-5. Accessed 2018-11-28
  • [31] Mason, D.: pydicom documentation (2008–).
  • [32] Oliphant, T.E.: Guide to NumPy. Trelgol Publishing (2006–)
  • [33] Negishi, R., Sekiguchi, K., Totsuka, Y., Uchida, M.: Determining Parastichy Numbers Using Discrete Fourier Transforms. Forma (2017). doi:10.5047/forma.2017.003
  • [34] Liew, S.F., Noh, H., Trevino, J., Negro, L.D., Cao, H.: Localized photonic band edge modes and orbital angular momenta of light in a golden-angle spiral. Optics Express 19(24), 23631–23642 (2011). doi:10.1364/OE.19.023631. Accessed 2018-06-28
  • [35] Eberly, D.: Fitting 3d Data with a Helix, 4 (2008)


Figure 1: \csentenceCherry tree trunk and 3D reconstruction. MRI produces a voxel-based image of the specimen. Using 3DSlicer, an open-source software tool, we created a 3D reconstruction. Two types of traces (a branching, or ”V”, trace and a single trace) are shown as physical cross-sections made by cutting the trunk and as slices from the scan.
Figure 2: \csentenceAnatomy of sweet cherry wood. A labeled slice from the MRI scan of sweet cherry, colored by intensity value. The pith is the hollow center of the trunk. Both the branch and the epicormic trace are high intensity regions.
Figure 3: \csentenceTwo views of parastichy spirals from artificial data. A set of 40 nodes, generated by the golden angle and constant rise (). In the centric form, parastichies are spiral arms emanate from the origin and in the cylindrical form, parastichies are helices in a cylindrical lattice, called the Bravais lattice. In the Bravais lattice, the and rays are the same, allowing the parastichies to “wrap around” the circumference of the plant.
Figure 4: \csentenceAlgorithm to determine the approximate tree radius in each slice. (A) To estimate the radius of the trunk in each slice we first apply the Sobel Edge Detection algorithm to reveal areas of abrupt change in pixel intensity. (B) Then we threshold the image with equal to the 99th percentile of the values given by edge detection. (C) Finally we compute the median distance from the pith location to all of the pixels in the binary image to estimate the radius. (D) To estimate the overall radius for the trunk, we compute the median of all radii for each slice. The dotted line represents the median
Figure 5: \csentencePolar coordinate conversion for selected slices. After finding the center and radius of each pixel slice we convert to polar coordinates, producing polar slices . The radius spans from 0 to the estimated radius .
Figure 6: \csentenceIsolating phyllotactic features in the boundary image. (A) The boundary image, created by assembling the means along each radial direction in the polar slices. Rows correspond to polar coordinates and columns slices. (B) The thresholded boundary image, created using Otsu’s method. (C) Cylindrical Bravais lattice of nodes, with 2- and 3-parastichies shown. The blob corresponding to node 26 was added by choosing a less restrictive threshold for that region.
Figure 7: \csentenceDistributions of phyllotactic parameters. (A) The distribution of local divergence angles compared to the golden angle () and the mean (). The intervals given by Adler’s Theorem (creftype 3) for visible and opposed parastichy pairs are shown below the distribution. Assuming the mean of the distribution is the true divergence angle, (1,2), (2,3), (3,5) and (5,8), are visible and opposed. The interval [, ], corresponding to the pair (1,2), is not shown. (B) The distribution of rise values, showing the range of internode lengths.
Figure 8: \csentenceVisualizing parastichies in 3D. Families of parastichies using helical interpolation on a cylinder [35]. (A) The blobs found in Fig. 6(B) map to wedge shaped regions in the MRI scan. (B) A 3D visualization of parastichy families that are a part of visible and opposed parastichy pairs. (C) Parastichy families shown in the Bravais lattice, with the same color scheme as part (B).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description