A fibre tracking algorithm for volumetric microstructural data - application to tendons

Tendons are crucial connective tissues made almost entirely of bundles of long, near-parallel collagen ﬁb-rils, and are vitally important to skeletal stability and mechanical function. Tendon structure is typically quantiﬁed in 2D, whereas, in this work, we have used serial block face-scanning electron microscopy to image tendons in 3D. We present a custom ﬁbre tracking algorithm (FTA), with which we have characterised the 3D microstructure of tendon. Currently available tools for ﬁbre tracing were unsuitable for tracking large numbers of ﬁbrils and handling imaging artefacts associated with EM. We have tracked ﬁbrils through a representative tendon volume and measured their relative length, diameter, orientation, chirality, tortuosity and volume fraction, which are just some of the measurements it is possible to make with the FTA depending on the research question. This algorithm has been developed in a general way and can be applied to a range of biological research questions relating to tendon structure-function relationships, on topics such as ageing, disease, development and injury. The FTA is also applicable to other ﬁbrous biological materials, as well as engineered materials and textiles; it is written using Python and is freely available to download.


Collagen fibrils in the extracellular matrix
The extracellular matrix (ECM) provides the scaffold to give structure and strength to biological tissue.ECM accounts for 70% of aligned network of fibrils in the skin gives rise to more isotropic stress tolerance [3,4] .
Altering the architecture of fibrous scaffolds can have a significant impact on the behaviour of cells and the ability of a tissue to function correctly.Changes to the matrix organisation affect cell migration [5] , stem cell fate [6] , cell morphology [7] and even matrix production [8] .Additionally, in some cases, cells are able to sense mechanical force or cell-matrix geometry, which is transduced into biochemical signals that result in cell responses [9] .Furthermore, because the ECM is often disrupted in cases of disease or injury in both heritable diseases like Ehlers-Danlos syndrome [10] and acquired diseases such as fibrosis or cancer [11] , measuring the arrangement of fibres within tissue samples can be used for diagnostic purposes and for monitoring the response of tissues to the appropriate clinical therapies.Therefore, characterisation of native matrix scaffolds or engineered fibrous materials from imaging data is crucial across a wide range of fields of interest within biomedical research disciplines and in understanding collagen related disorders.

Tendons are a model fibrous tissue
Tendons connect muscles to bones across joints and thus enable forces generated by muscle contractions to initiate skeletal movement.The stability and mechanical function of animals depend on the healthy functioning of tendons and other connective tissues which carry tensile and compressive loads.Tendons are composed of bundles, or fascicles, of long, predominantly type-I collagen fibres surrounded by an extrafibrillar matrix containing mainly water and proteoglycans.In biomechanical terms, the tendon is often regarded as a fibre reinforced composite [12,13] .Tendons are largely acellular, the volume fraction of cells compared to ECM is small, measured at 2% in adult rats [14] .Additionally, the elastic modulus of cells is orders of magnitude lower than that of tendon fibrils.The equilibrium Youngs modulus of chondrocytes (cells found in cartilage) is around 0.2 kPa, which is smaller than the longitudinal Young's modulus of a typical tendon fibril by around six orders of magnitude.Tenocytes (cells found in tendon) may be stiffer, but not by an amount that could make them mechanically significant in comparison to matrix elements [13] .
Fibrils exhibit a regular, wavy, longitudinal morphology called crimp, thought to undulate in either a two-dimensional (2D) planar or helical fashion [15][16][17] over an approximately 100 μm period for adult mouse tendon [18] .Their diameter remains constant along most of their length, with some exceptions [19] , but there is a distribution of diameters in any given tendon cross-section.Recent work has shown that the diameter distribution can be described well by a trimodal Gaussian distribution [19] .
It has been noted previously that gross tendon mechanical properties are highly dependent on the sub-fascicle microstructure and fibril sizes.An increase of 50% in fibril diameter has previously been found to account for a five-fold increase in ultimate tensile strength (the maximum stress before failure) and Young's modulus in human fibroblast tendon constructs [20] .Collagen fibrils in tendon are structurally continuous through large longitudinal distances ( ∼ 10 −2 m) [21] , and can be as small as 50 nm in diameter, and as large as 500 nm [19] , giving them an aspect ratio of 10 5 -10 6 .
Tendon is a useful model tissue for the study of a range of diseases and heritable collagen disorders.For example, Sun et al. [22] have created a mouse model with a collagen V mutation similar to that associated with classical Ehlers-Danlos Syndrome.The fibrils in these mouse tendons were observed to be less regularly packed, with more space separating them than in the corresponding controls in the study [22] , and other groups have noted altered fibril morphology and diameter distributions in similar mod-els compared to wild-type tendons [23] .Mice with osteogenesis imperfecta type phenotypes often have tendons of a smaller size and show irregular fibril alignment [24] , so the accurate quantification of three-dimensional (3D) tendon microstructure is important in this field.

Challenges in imaging collagen fibrils and fibrous tissue
Imaging fibrous tissue can be technically challenging.Generally, the success of an imaging method depends upon the information desired in the experiment, the fibre sizes, and the fibre density.For collagen fibrils, the high aspect ratio provides a challenge in imaging since the fibril width renders optical microscopy unusable (such techniques are bound by a diffraction limit of approximately 10 −7 m).Super resolution techniques can break this limit, but they rely on fluorescence and antibody staining [25] , which can only be used over a limited z-depth, meaning these techniques are unsuitable for imaging long tendon fibrils.
Imaging of neurons and, to a lesser extent, fibrous tissues, has been achieved using diffusion-tensor magnetic resonance imaging, but this approach is naturally limited to structures that can be resolved using magnetic resonance [26][27][28] .It is also possible to use polarised light microscopy to measure fibre alignment, but this technique is also diffraction limited and further requires that the materials imaged are birefringent, examples of which include collagen and nerve fibres [29][30][31] .X-ray computed tomography (XCT) can be used to visualise 3D structures, and is non destructive; however, the resolution is insufficient to see collagen fibrils in cross-section, though it has been used to visualise bone, tendon and ligament on a larger scale [32][33][34][35] .
For fibrils in tendon, scanning electron microscopy (SEM) and transmission electron microscopy (TEM) are used due to the diffraction limit of light microscopy, and yield exceptional results.These techniques are commonly used to produce 2D planar images to give a representative slice perpendicular to the tendon axis [19,36,37] ; however, tendon function depends not only on fibril diameter distributions but also on the number of fibrils, their 3D organisation, volume fraction, crimp structure and length.Using a single 2D micrograph to gather only a fibril diameter distribution ignores many of the intricacies of global fibril architecture and loses information about lateral variations in fibril diameter and longitudinal variations.
Serial block face-SEM (SBF-SEM) is an extension of conventional SEM, and uses a microtome to remove thin slices of a fixed and embedded sample [38] .The top of the block face is imaged before the sample is cut at a specified depth by the microtome, revealing a new block face, which is again imaged and the process repeated for as great a depth as desired [39] .SBF-SEM makes it possible to track tendon longitudinally whilst capturing individual fibrils in cross-section.Previous work has documented fibril orientation and mean fibril lengths using SBF-SEM on small volumes [38,40] and SBF-SEM is more technically demanding and time consuming than conventional SEM and TEM, but can illuminate the 3D structure of tendon in a way that is not possible with these techniques.

Automatic fibre tracking
Identifying fibres from images is sometimes conducted by manual segmentation.This is where pixels are assigned to objects individually by tracing over them by hand in the images.This is done in both in 2D and 3D [41] , but the 3D case is much more timeconsuming.Moving beyond manual characterisation, there are a number of open-source tools that can be used to measure various fibre properties in 2D.For example, the commonly used FIJI package has various plugins which are useful for determining fibre orientations.These include FibrilTool and OrientationJ [42,43] .
Elsewhere, CytoSpectre is a tool that can evaluate orientation and size distributions of structures in microscopy images [44] .Unfortunately, moving from analysis of 2D structures to 3D structures using such tools can be problematic.Fibres oriented approximately perpendicular to the image plane may be difficult to detect, and information about their longitudinal morphology can be lost.Additionally, the calculation of fibre diameter is not always possible directly from 2D images as the apparent cross-section can be misinterpreted depending on the angle of orientation with respect to the image plane.
The complexity introduced by 3D image analysis presents a challenge and limited options exist for 3D quantification of fibre geometries.Fabric tensors can be used to measure fibre orientations and directions of anisotropy in porous materials, but are not suitable for fibrous tissue [45][46][47] , structure tensors can be used to measure mean orientations in 3D, but cannot be used to track individual structures [48] , and the Fast Fourier transform, which is commonly used in 2D, has been adapted to 3D by comparing Fourier transforms with a filter bank to calculate the orientation of collagen in second-harmonic generation images [49,50] .This analysis only provides bulk orientation information, however, and is therefore unsuitable for resolving individual collagen fibrils.
The development of a vector summation technique [51,52] has been successful in producing 3D quantification of fibre orientation and individual fibre analysis.The fibre orientations, however, are accurate only to 4-5 • , on average, and the technique cannot measure fibre diameters or account for branching.A fibre centreline tracing method was developed to segment fibres and measure their diameter, length and orientation [53,54] .This was then further developed into the curvelet transform fibre extraction (CT-FIRE) application [55] .This application has limitations, however, including difficulty in identifying highly curved fibres in dense networks or those with variable thickness and it is computationally expensive [56][57][58] .Furthermore, the current, publicly available version of the application is limited to 2D, and is designed for networks of meshed fibres with many junctions, which is not particularly useful for long, mostly aligned parallel fibres such as those in tendons.
In the fields of materials science and textiles, there have been attempts to trace long parallel fibres through reinforced composites in order to find their orientations.These methods track the centre points of the fibres in successive planes [59] , sometimes with a Kalman filter [60] .Since fibres in engineered materials are usually of a known diameter (and are thus not measured in the analysis), and are very straight and parallel, these approaches are not suitable for direct use on biological tissues, which have a more disordered, noisier microstructure.More recent work [61] has quantified orientation and diameters for fibrous materials in a voxel-wise fashion, but struggled with high density fibre populations, meaning they were unable to segment collagen fibrils in tendon, and instead identified sparse elastin fibrils [61] .Thus, while there are tools available for automated 3D measurement of fibrils in tendon, current techniques are limited in their accuracy or the extent of their capabilities.

Our approach
We have formulated a method for automated fibril tracking to quantify microstructural details in SBF-SEM images of tendon volumes using a substantial amount of imaging data that captures a wide field of view in three dimensions.Fibrils are tracked through the volume, producing a database containing centroid and geometric information for each fibril in each plane.There are many possible measurements which can be calculated using this database, for a range of research questions.Our algorithm and its out-puts can be applied within tendon research to questions of tissue mechanics and used to identify morphological differences between tissues of different species, ages or disease states.To illustrate some of the capabilities of the fibril tracking algorithm (FTA), we present measurements of individual fibril diameters, relative lengths, alignments, tortuosity and chirality.The algorithm is developed in a general way and is applicable to images of any material with aligned fibres (such as ligaments, blood vessels, engineered materials and textiles) captured using any imaging modality.

Biological data acquisition
This study is based on archived SBF-SEM data collected as part of previous research [19] and therefore no new animals were sacrificed.The data was collected at the University of Manchester in compliance with institutional ethical protocols.The Achilles tendon of a wild-type C57BL/6 mouse was dissected, embedded in resin and stained with osmium tetroxide in line with standard procedures for imaging with SBF-SEM.The process is outlined in full by Starborg et al. [38] .The care and use of all mice was in accordance with UK Home Office regulations, and the UK Animals (Scientific Procedures) Act of 1986 under the Home Office Licence (#70/8858 or I045CA465).

Image processing
The required input for the FTA is a stack of accurately segmented images, so it is necessary to do some pre-processing to handle the raw greyscale biological images.Collagen fibril sections in tendon, when captured by SBF-SEM, are a roughly consistent grey value due to their near-uniform molecular density.They are closely packed, mostly circular, and of variable diameter.Collagen fibrils make up the vast majority of tendon tissue by volume, but cellular matter and other matrix proteins are also present, and are often seen in micrographs ( Fig. 1 A).This provides us with a specific set of challenges for image segmentation.Of course, in a different tissue or material, these challenges may be different, but here we provide a flexible image segmentation pipeline using Cell-Profiler [62] .This allows efficient processing of large, consecutive batches of images.
The pipeline is summarised in Table 1 , and the associated Cell-Profiler script is available in the Research Data we have made available ( Section 5 ).The micrographs are processed using this pipeline so that the fibrils are separated and much of the cellular matter is removed, as our aim is to track fibrils only.A typical processed image is shown in Fig. 1 B. It is necessary to inspect the volume to check for any broken planes.For each broken plane, the image number is noted and used as an input into the algorithm, so they are skipped over when tracing fibrils in order to avoid erroneous results.These damaged planes are an artefact of SBF-SEM, which we discuss further in Section 4 .

Overview of tracking method
Our algorithm works by measuring properties of the objects in plane p (these objects are denoted by an index i , ranging from 0 to n p − 1 , where n p is the number of objects in plane p) and comparing them to the objects in plane p + 1 (denoted by index j, ranging from 0 to n p+1 − 1 ).A simplified representation of this process is given in Fig. 1 C.A set of criteria based on the centroid location, area and minimum Feret diameter (MFD) of each object is used to match the objects through successive planes.These properties of  1) ) are generated for a pair of planes.For this simplified example of three fibrils in planes p and p + 1 , the 3 × 3 matrix of the values of ξ (p, i, j) would contain smaller values along the leading diagonal than the off-diagonal elements, as these are the correct (i, j) pairings.D : The red fibril segment i highlighted in plane p is compared to all the segments j in plane p + 1 in the white square, including those touching the border, and each value of ξ (p, i, j) is calculated and added to the error function matrix.This is repeated for every fibril segment i in plane p to evaluate the best matches for the segments in planes p and p + 1 , then repeated for all planes in the volume.The white square has side length equal to 10 times the MFD of the fibril being tracked (single, red fibril, left), and is centred on its centroid.E : Measuring fibril chirality.Starting from the segment at the top of the image stack, the trajectory in the x − y plane of a helical fibril is plotted.Since a right-handed helix will twist in a clockwise fashion, the cross product can be used to determine the chirality of the fibril.The position vectors of three consecutive fibril segments (u 1 , u 2 , u 3 ) are used to calculate the cross product in Eq. ( 13) .The only non-zero entry of the resulting vector is always its third entry as the position vectors are co-planar in the x − y plane.For a clockwise path between three adjacent fibril segments, this value will always be negative and for an anticlockwise path, it will be positive.The sign of this cross product is calculated along the whole length of the fibril to quantify its chirality (see Eq. ( 13) ).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)the fibril segments were measured using Scikit-image [63] , with the exception of the MFD, which was calculated using the convex hull of each object and a custom script (included as part of the FTA -see the Research Data statement for details on how to access the FTA).

Error minimisation method
To join fibril segments together, we use a metric we call the error function, which is a measure of the likelihood that two fibril segments i and j in planes p and p + 1 , respectively, are part of the

Gaussian filter
Smooth and remove imaging artefacts σ = 1 pixel.This is the standard deviation of the kernel to be used for blurring.A larger value of sigma induces more blurring.Thresholding Separate objects of interest from the background Adaptive thresholding was used, with a window size of 40 pixels.This was to account for uneven illumination in the image.The thresholding algorithm tests for a value where the cross-entropy between the foreground and the background is minimised.Remove objects that are too large or too small

Remove cells
Cells sometimes get split into many pieces in the thresholding step, so, by filtering for objects of width larger than 10 nm and smaller than 400 nm, the majority of the cell debris is removed.The fibrils are still mostly touching at this point, so no fibrils are removed at this stage.

Opening
Smooth the edges of the fibrils Morphological opening allows the spaces between the fibrils to be revealed, and removes small bridging artefacts and smooths the fibril edges after the thresholding step.A disk matrix of size 3 pixels was used.

Watershedding
Separate touching fibrils Local minima (the centroids of the fibrils) are calculated, and then the image is treated topographically.The image is 'flooded', from the local minima, and the fibril boundaries are introduced at the points where the pools touch [63] .

Remove border objects
Remove fibrils which are not fully in the field of view We do not wish to include fibrils which are bisected by the field of view to ensure that only whole fibril segments are tracked.
same continuous fibril.The error function ξ is defined as a weighted sum of three contributions, ξ c , ξ d , and ξ a , which quantify the differences in the centroid (c) , MFD (d) and area (a ) , respectively, of segments i and j.The values of ξ c , ξ d , and ξ a are weighted by the constants α, β, and γ , respectively, and ξ is normalised by their sum.The default values of the parameters are ) and these values were used to produce the results in this paper.The selection of these values is discussed in Section 3.2 .The error function has been designed so that the most realistic candidates for matching will yield an error ξ (p, i, j) < 1 .
We use a coordinate system which defines the x -and y -axes as being in the plane of each cross-sectional image and the z-axis as following the axis of the image stack, which is broadly aligned with the tendon axis.

Centroid error: ξ c
The centroid is the most important component of ξ (p, i, j) as two consecutive cross-sections of the same fibril will be at a similar location in x and y .To calculate ξ c , the Euclidian distance between the centroids of the i th object in plane p and the jth object in plane p + 1 is calculated and normalised by a characteristic length scale for the data.Thus, ξ c is given by where .denotes the Euclidian distance (L2 norm) between the centroids of fibril i in plane p, (x p i , y p i ) , and fibril j in plane p + 1 , (x p+1 j , y p+1 j ) , and L i is a chosen length scale.This length scale should be selected by the user to be the maximum likely in-plane distance between two consecutive cross-sections of the same fibre.For the SBF-SEM tendon data we analysed in this paper, we selected this value to be the MFD of object i , D p i .The value of L i can be chosen differently for other data sets by manual inspection before running the FTA.We have chosen to vary the value of L i depending on the object being considered, but the user could choose to use the same value for all objects if required.Careful selection of L i ensures that, for realistic object pairings, the value of ξ c (p, i, j) is less than or equal to 1.

Minimum Feret diameter and area errors: ξ d , ξ a
The MFD and area are both included in the error function ( Eq. ( 1) ) to ensure that fibril matching accounts for the size and shape of each segment as well as its location.This is necessary because of cases when, due to fibril undulation and artefacts caused by sample sectioning, the segment whose centroid is nearest to the target segment in the previous slice is actually from a different fibril.It is clear by eye that such cross-sections should not be matched, but an algorithm based on centroid location alone would erroneously do so.
The errors ξ d (p, i, j) , ξ a (p, i, j) , which correspond to the MFD and area, respectively, are calculated in the same way, by quantifying the differences in these quantities for two objects i and j.They are expressed as and where D p i and A p i are the MFD and area of fibril i in plane p, respectively.These errors thus represent fractional changes in MFD and area.The errors ξ f (p, i, j) and ξ a (p, i, j) will always be positive and for a good match of an i , j pair, they will be small.

Details of the tracking algorithm
In the first plane ( p = 0 ), an n 0 × n 1 array (where n 0 and n 1 are the numbers of objects present in planes 0 and 1, respectively) of ξ values is created.To improve computational efficiency, only objects j which lie within a square with side length equal to 10 times the size (MFD) of the object i and centred on its centroid are used to populate the array.This is represented visually in Fig. 1 D. The ξ values are then sorted in ascending order, and the corresponding i , j pairs listed in order of increasing ξ .An iterative process of assigning pairs is then conducted as follows: 1. the pair with the lowest value of ξ , (i 0 , j 0 ) , is assigned as a match, 2. the object j 0 is appended to the entry in the fibril record which contains object i 0 , 3. every other (i, j) pair where i = i 0 or j = j 0 is removed from the list.
This process is then repeated, each time removing unwanted potential matches from the sorted list at each step.Not all segments in plane p will have a match in plane p + 1 .Some fibrils leave or enter the field of view, or merge to a neighbouring fibril, or stop entirely.After a number of iterations, all of the true matches will have been made between plane p and plane p + 1 , and the remaining values of ξ in the list will be relatively large (corresponding to unfit matches).We set a threshold value of ξ ( ξ lim ) which corresponds to the point at which we consider all further matches to be inaccurate.Here that value is set at ξ = 1 (we discuss this choice in Section 3.4 ).
At this point, we establish a record of our joined segments, we mark any unmatched segments in plane 0 as discontinued, and add any unmatched segments in plane 1 to the end of the fibril record, denoting them as potential starts of new fibrils to be traced in subsequent planes using the same process.This process is repeated through successive planes until all the good matches ( ξ < ξ lim ) have been made in all planes.

Selecting fibrils of interest
As mentioned previously, fibrils are regarded to be structurally continuous throughout tendons, though in the sample imaged by SBF-SEM, we may not see fibrils that are continuous throughout the whole volume.This is because the field of view travels in a fixed z-direction (which is perpendicular to the plane of sectioning), and the fibrils themselves are not straight, as they are crimped.Additionally, the sample will likely not have been aligned well enough for the longitudinal axis of the tendon to perfectly co-align with the z-axis at the image stack.As a consequence, many of the fibrils seen in the first plane leave the field of view in a subsequent plane.Therefore, our algorithm produces a fibril record which catalogues many joined segments of various lengths.Depending on the application, at this point, a subset of the fibril record can be extracted for analysis.As an example, in this paper, we select fibrils which appear through a continuous distance of at least 10 μm in the z-direction, and discard the remainder.Once the population of interest has been selected, it can be measured.The raw output of the algorithm is a database, with entries corresponding to individual fibrils, which records their centroid co-ordinates, MFD, cross-sectional area and eccentricity (measured by fitting an ellipse) of each fibril in each plane.

Measuring the tracked fibril population
Below we provide some examples of quantities we have chosen to measure.We report the fibril relative length, MFD (as a mean value along each fibril's length), alignment with respect to the fascicle in which it is embedded, chirality and tortuosity.Using the fibril record generated at the end of the procedure described in Section 2.5 , it is also possible to measure the variation of MFD or area along the length of the fibres, as well as their eccentricity using an elliptical fit, for example.With further computations using the centroid values of the fibril population, a study into nearest neighbours or fibril separation distance could also be conducted, both in-plane and variably along the z-direction.As it is adaptable to the research question being asked, the FTA is able to generate a wide range of metrics of the tendon microstructure.

Fibril length and tortuosity calculations
Across the tracked fibril population, we calculate the relative length λ i of fibril i , which is a measure of how long an individual fibril is compared to the fascicle in which it is embedded.This is an important parameter when predicting the mechanical behaviour of a tendon [64][65][66] .To calculate this, we calculate a parameter we call the equivalent fascicle length .The equivalent fascicle length is calculated from the mean movement of the fibrils in each plane as follows: where ( x p , ȳ p ) is the average in-plane centroid location of all the tracked fibrils in plane p, and z p is the z-coordinate of plane p.We sum between plane p i 1 and p i 2 , the first and last planes in which fibril i appears.Using this quantity we calculate λ i as where l i is the tracked length of fibril i , which is calculated by summing the Euclidean distances between consecutive centroids along the fibril's arc length: where (x i p , y i p , z i p ) is the centroid location of fibril i in plane p.The tortuosity τ i of fibril i was calculated as where, l SL i is the straight-line distance between the centroid of the first and last fibril segment.A greater difference between these two lengths increases τ i and indicates a more tortuous fibril.

Fibril alignment calculations
Fibril alignment is calculated by determining the angle between the fibril and the fascicle in each plane.We define the unit vector ˆ v i p , which gives the direction of a particular fibril i in a particular plane p, as where r i p is the position vector of fibril i in plane p.We similarly define another unit vector ˆ v f p , which describes the direction of the fascicle in plane p as where r p is the mean position vector of all the traced fibrils in plane p and r p .The alignment θ i of fibril i with respect to the fascicle is calculated as the mean angle between fibril i and the fascicle over all the planes in which fibril i is present: 2.6.3.Fibril chirality calculations Fibrils have been observed to follow helical paths within tendon [67] .To quantify the helical nature of fibrils, in this work, the fibril paths in the x − y plane were used, collapsing the z-direction.If the 2D fibril path spirals in a clockwise fashion, the helix is righthanded, and it is left-handed if the fibril path spirals in an anticlockwise fashion.Using a Savitzky-Golay filter to smooth the 2D fibril path, the following cross product was calculated, where C i k is a quantity that is related to the local chirality of fibril i at point k , for 1 < k < N i , where N i is the number of points that make up the smoothed path of fibril i and u i k is the position vector of point k , which will have 0 as its third entry since all points have been collapsed onto a single plane.The sign of C i k indicates the local direction of curvature: when C i k < 0 , the path is curving in a clockwise fashion locally, and when C i k > 0 , it is curving in an anticlockwise direction.This is illustrated in Fig. 1 E. The number of positive and negative values of C i k are counted for each fibril.A fibril which is helical will amass an imbalance of positive or negative values, whereas a straight or non-helical curved fibril will have roughly equal numbers of positive and negative values of C i k .A chirality parameter H i , which quantifies how helical the fibril is, was calculated for each fibril as follows: where sgn is the sign function.The fibril paths in the x − y plane were inspected manually to discern a threshold value of H i to distinguish between helical and non-helical fibrils.It was observed that fibrils whose chirality parameter satisfies H i < −0 .15 were right-handed, and those whose chirality parameter satisfies H i > +0 .15 were identified as left-handed; the remainder (with −0 .15 ≤ H i ≤ 0 .15 ) were regarded as non-helical.

Results
The FTA provides a novel method for automatically segmenting fibrils in volumetric SBF-SEM image data (this approach could also be applied to other 3D imaging methods such as XCT).Some example results are shown in Fig. 2 .The aim was to trace tendon fibrils in 3D and gather corresponding morphological data from the reconstructed volume.A colourised volume displaying the full 3D reconstruction of a tendon sample is shown in Fig. 2 A and an animation moving through each plane is shown in video 1.This demonstrates the ability of the algorithm to trace large numbers of fibrils through a considerable distance in z, which was one of the primary aims of this work.

Observed microstructural measurements
Selected microstructural information derived from the reconstructed volume in Fig. 2 A is plotted in Fig. 2 B-D.Fibril length relative to the fascicle length is plotted in Fig. 2 C).We have fitted a normal distribution to this data, which gives a mean relative fibril length of 1.034, illustrating that the fibrils in this volume are, on average, 3.4% longer than the fascicle in which they are embedded.Another important measurement is that of the diameter distribution ( Fig. 2 C), which is important when studying healthy tendon development, ageing and disease as mentioned in Section 1.1 .To this data, we fitted a bimodal probability density function (PDF), which consists of the weighted sum of two normal distributions: where N is the normal PDF, w (which is confined to the range 0 < w < 1 ) is the weighting factor, μ 1 and μ 2 are the means of the two normal distributions and σ 1 and σ 2 are the associated standard deviations.This composite distribution fits the data reasonably well and highlights two populations of collagen fibrils.The first population has an MFD of 138 nm and a standard deviation of 44.1 nm and the second has a mean of 227 nm and a standard deviation of 17.1 nm, the weighting factor, w = 0 .7 .Our data show that the fibrils are also closely aligned with the longitudinal axis of the fascicle.A histogram showing the frequency of each angle of alignment, as defined in Eq. (11) , is plotted in Fig. 2 D. We have fitted a log-normal distribution to the data as θ is strictly positive, which gives a mean alignment of 6 • relative to the fascicle.
Using the method described in Section 2.6.3 , the chirality of each fibril was calculated using Eq. ( 13) .The spread of H i values is given in Fig. 2 E. These data were fitted with a skew normal distribution with location parameter μ = 0 .041 , scale σ = 0 .147 and shape a = −1 .28 .We found that 883 of the 4075 fibrils (22%) in this image volume were helical, of which 763 (86%) are righthanded and 120 (14%) are left-handed.Safa et al. found tendon fibrils to have 62% right handed helices [69] , though they used a much smaller sample size of 21 fibrils.
The calculated fibril tortuosity is given in Fig. 2 F. These results display greater tortuosity than Safa et al. found [67] .They reported that τ is typically less than 1%, whereas, here, we found the distribution of τ to be centred around 5%.To these data, we have fitted a log-normal distribution, as the values of tortuosity are all positive ( Fig. 2 F).We did not find any notable correlation of tortuosity with fibril diameter for helical fibrils or non-helical fibrils, which does not agree with the findings of Safa et al., who reported a negative correlation of tortuosity with increasing fibril diameter in helical fibrils.These findings may arise from the much larger sample size in this work, the difference in how the chirality was determined and calculated, and the fact that this image volume is from an energy-storing (Achilles) and not a positional (tail) tendon, as in Safa et al. [67] .

Accuracy of the algorithm
The accuracy of the tracking algorithm is visible in video 1, where each fibril track is simultaneously colourised.Segments which belong to tracks that are too short, and are thus not included in the final population are white in this video.Fig. 3 A illustrates the accuracy of the algorithm in the central plane of the volume.As can be seen, the vast majority of fibrils are captured (shown in red), whereas only a small fraction of the fibrils and some other objects (shown in grey) are discarded.Some of these objects are cellular debris and are excluded on the basis of their morphology (since they are far from circular).
To test the robustness of the algorithm, we compared the MFD of the fibrils which are tracked to the MFDs of all fibrils in a few sample planes which were hand-selected to ensure that no cellular matter was present.This was done to ascertain whether the tracked fibrils are representative of the whole fibril population.The tracked fibrils, as mentioned in Section 2.5.1 , are fibrils that span more than 10 μm in the z-direction, and all other tracks are excluded.Fig. 3 B shows the MFD distributions plotted together.As seen in the plot, the FTA tracks large fibrils with high accuracy, picking up more than 90% of those with an MFD above 150 nm, but is less successful at tracking small fibrils, finding only 60% of them with MFD under 120 nm.One reason for this is that smaller fibrils can get lost at the image segmentation step, either by being grouped with a nearby large fibril, or because they have insufficient contrast to be detected in some planes.This can mean that even if the fibril is traced, the track is disjointed and discarded on the basis of not meeting the 10 μm z-depth condition.A multimodal distribution of fibril diameters is seen here, in agreement with previous work [19,68] .The red line is a bimodal normal probability density function fitted to the data.D : The orientation of the fibrils relative to the fascicle.The red line is a log-normal probability density function fitted to this data.E : The chirality of the fibrils, calculated using Eq. ( 13) .The dashed, black line is a skew normal probability density function fitted to the data, with location, scale and shape parameters μ, σ , a .Fibrils which had a chirality value greater than 0.15 or less than −0 .15 were observed to be helical, by inspection.F : Fibril tortuosity, calculated using Eq. ( 8) .The fit given is a log-normal distribution.The parameters in the legend for D and F are the log-normal distribution parameters, which represent the mean and standard deviation of the variable's natural logarithm, not the mean and standard deviation of the variable itself.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Depending on the research question being asked, this may not be too important.60% is still the majority of fibrils in this population, and they are the least mechanically relevant.If the proportion of lost, small fibrils is consistent across different samples, it will still be possible to compare them bearing this limitation in mind.

Handling of imaging artefacts
One of the issues raised by SBF-SEM specifically is artefacts in sectioning which can cause the whole sample to slip in a particular direction.An example slipping artefact is displayed in Fig. 3 C.A sample slip such as this renders this plane unusable.In a typi- The MFD of the tracked fibrils (blue) and the whole fibril population in a randomly selected 2D plane (grey).The dashed, black line represents the proportion of fibrils of each size that were successfully tracked, which ranges from approximately 60% to 100%.The FTA successfully tracks fibrils with an MFD of 150-300 nm, but is less successful in tracking fibrils with an MFD of less than 150 nm.C : A typical broken plane.This is seen where the sectioning is unsuccessful.In a volume of 700 planes, there are typically fewer than 10 broken planes like this.It is important that the FTA correctly handles these artefacts, otherwise erroneous matches or discontinuities in the fibril tracks will occur.D : The response of the FTA to broken planes.The arrows indicate the locations of broken planes within the volume, and the vertical axis indicates the number of fibril termini which are not close to the edge of the image (they are more than 5% of the width of the image from the border to exclude fibrils that leave the field of view from the calculation).There are typically fewer than 10 termini per plane, which, in principle, could either be true fibril ends or an artefact of tracking.Broken planes do not tend to lead to a large increase in termini, indicating that the FTA is handling them successfully.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)cal SBF-SEM volume, such as the ones used in this study, typically, fewer than 1 in 100 planes are affected in this way, but they can cause errors in the FTA.Here, we have dealt with this via an input to the algorithm which lets the program know which planes to skip over.It is thus necessary to check through the volume and manually input any unusable planes to the FTA.To test the impact these broken planes have on the overall tracking, Fig. 3 D shows how many fibril ends there are in each plane.In the first and last plane, we see the most, since the fibrils in the volume are largely continuous through all 700 planes.At each plane, we see a certain number of fibrils that terminate there.Some of these will be erroneous, but they are mostly fibrils either entering or leaving the field of view.If the algorithm was poor at handling broken planes, we would see a large increase in fibril ends after a damaged plane, as the algorithm would be unable to find appropriate matches.This is generally not the case.In Fig. 3 D, the locations of damaged planes are represented by arrows, and it is evident that, with the exception of at the black arrow, where several planes were damaged relatively close to each other, there is not a noticeable increase in fibril ends before or after them.Even at the black arrow, only 22 fibrils were lost, which is a very small number compared to the total number of fibrils that were successfully tracked through this slice.This shows that the algorithm is robust enough to handle this kind of artefact.

The error function
When matching consecutive fibril segments, the error threshold determines which matches are accepted.If we set this threshold too high, inaccurate matches are made, whereas if we set it too low, not all tracks are properly captured.The value of the error threshold used in this work is based on the maximum likely values of ξ c (p, i, j) , ξ d (p, i, j) , ξ a (p, i, j) being equal to 1. Since we have set the constants (α, β, γ ) = (1 , 1 , 1) , the default error threshold is also ξ = 1 by Eq. ( 1) .The effect of using this threshold value is displayed in Fig. 4 A. Fibril segments from a randomly selected subset of those that appear in a randomly selected plane are each represented by a coloured line, and their MFDs are given in the legend.For each fibril segment, the potential candidates to match to in the next plane are ranked.The numbers on the horizontal axis represent the ranking of each candidate and the corresponding error is given on the vertical axis.The error threshold which we have set is represented by the dashed line at ξ = 1 .As can be seen, all but one, small segment (MFD = 61 nm) has several suit- able matches, whereas for the small segment (the darker blue line), even the candidate that is ranked first produces a high error, so no matching will occur.This shows that the FTA is good at finding matches for all but very small fibrils, which is discussed below.
The error minimisation equation used in the FTA ( Eq. ( 1) ) is designed to be flexible regarding the importance of centroid location, MFD and area when matching consecutive segments.We designed the three errors ξ c (p, i, j) , ξ d (p, i, j) , ξ a (p, i, j) to be nondimensional and the weights associated with these errors, α, β, γ , can be varied to ascertain whether greater accuracy can be achieved in the parameter space around the point (α, β, γ ) = (1 , 1 , 1) .The results of varying these parameters are plotted in Fig. 4 B. This heat map represents the fraction of failed matches between randomly selected planes.It is generated by setting α = 1 , whilst β and γ are varied from 0 to 3. Since the three constants are normalised by their sum ( Eq. ( 1) ), it is the relative values of α, β, γ that are important and not the values themselves, which is why α remains fixed.In the parameter space where α, β, γ > 0 , there is little variation on matching frequency, as shown by the narrow range of values (2-3%) in this plot.Nevertheless, to investigate the effect that varying these parameters has on matching, we have presented the results of tracking for three different points in the parameter space, for a random sample of 50 fibrils for ease of viewing ( Fig. 4 C(i)-C(iii)).The first image presented here (i) corresponds to our starting values, (α, β, γ ) = (1 , 1 , 1) , and shows continuous fibrils, demonstrating that the default parameter values capture fibril tracks with high accuracy.The second image (ii) shows what happens when we set (α, β, γ ) = (1 , 0 , 0) , so that only centroid information, and no geometrical information (MFD or area), is incorporated.This leads to errors, as matches are made which select the nearest objects and effectively skip over to a different fibril track, leading to both inaccurate continuous tracks and disjointed short ones.The third image (iii) corresponds to (α, β, γ ) = (1 , 1 .5 , 1 .8) , which represents the greatest number of matches per plane (the darkest region).This mapping is reasonably accurate most of the time, but still leads to errors, as this parameter combination favours the morphology of fibril segments in the next plane over their location, which leads to erroneous jumps between clearly distinct fibril tracks.This investigation shows the design of the algorithm at (α, β, γ ) = (1 , 1 , 1) is fit for purpose, and that it is robustly and accurately tracking fibrils.For other applications, there is the flexibility to vary these three parameters depending on the microstructure of the material being studied.For example, to track fibres that each have a different cross-sectional size and shape or that are obliquely aligned, it would be beneficial to increase the values of the morphology parameters β and γ , whereas for fibres with similar cross-sections, these values could be lowered.

Limitations of SBF-SEM
This study was completed with SBF-SEM, which is a common technique for imaging structures that are of a size below the diffraction limit of optical microscopy.To prepare samples for SBF-SEM, tissue is chemically fixed, dehydrated, embedded in resin and baked at a high temperature (example protocol available in Starborg et al. [38] ).This process can be damaging to tissue and causes shrinkage due to dehydration (measured at 19% in lung tissue [70] ).It has been shown that some kinds of dehydration and fixation can cause reorganisation of collagen structure, both affecting the D-period and alignment [71] .This means that there will be a discrepancy between recorded microstructural information and the true structure in vivo when using images acquired in this way; therefore, measurements using this technique should be considered as indicative of the in vivo structure, but not quantitatively precise.For 3D work measuring objects of the size of collagen fibrils, SBF-SEM is the best and most commonly used imaging technique currently available and tissue shrinkage affects all studies that use this imaging technology in the same way.Some cryoelectron microscopy (EM) techniques such as those of Keene et al. [72] are able to preserve microstructure with higher fidelity, but these can currently only produce 2D images.We expect measurements of interfibrillar spacing, in particular, to be affected as the extra-collagenous phase is more greatly impacted by the dehydration step [70] , though there are no existing techniques to quantify and validate precisely how shrinkage affects each phase individually.If such measurements are reported in the future, it may be possible to apply non-affine transformations to the images to expand the fibrils and interfibrillar space to reflect their true morphologies in vivo .We note that, if a sample preparation protocol is developed in the future which ensures uniform shrinkage across different phases of the material, then relative microstructural measurements, such as the relative fibril length distributions measured in this paper will be unaffected by shrinkage.Furthermore, in spite of the issues associated with sample shrinkage, since all samples studied using SBF-SEM are affected in the same way, comparisons between different samples are still meaningful.
We have shown that our algorithm robustly deals with the imaging artefacts caused by SBF-SEM ( Section 3.3 ) which occur because this technique is destructive.As this is the most common method for resolving the 3D structure of collagenous tissues, it was important for us to account for these artefacts when designing the algorithm so it is applicable to these contexts.The FTA, however, is agnostic to the method of imaging.In a study of parallel fibres captured using a non-destructive imaging technique without damaged planes, the algorithm will run without skipping over any planes if the user does not specify to.This applies to XCT data, which is commonly used for 3D imaging of materials with larger microstructural features than those seen in connective tissue.

Image segmentation
As with any computational work using biological images, the quality of analysis is dependent on that of the images and segmentation.If they are of insufficient contrast or resolution, or contain sectioning or staining artefacts, the subsequent analysis can be erroneous.In tendons, such issues can make it difficult to distinguish individual fibrils from one another, or can create uncertainty in their sizes if their edges are not sufficiently sharp.This can occur if the images are taken at too low magnification, or if the staining and processing of the tissue itself is unsuccessful.If the stain fails to penetrate the tissue properly, the images will have insufficient contrast.Consequently, it is important to ensure that sample preparation and imaging is undertaken with precision and skill by specialists, as we have done in this study.When tracking tendon fibrils, it is of interest to ensure that cells and nonfibrillar objects are excluded from the analysis.Here, we have removed cells in the segmentation step in the CellProfiler script (see Table 1 and the script provided as described in the Research Data statement).This is based on their size and shape compared to the fibrils.Some fragments of cells may still remain at this stage, but are necessarily discarded in the tracking process, as they are not aligned or present through enough consecutive planes to meet the criterion we have specified that traced objects must span 10 μm in the z-direction.To improve efficiency in the future, a mask could be applied to remove the cells from the volume before segmentation, to avoid fragments of cells remaining in the image volume.We have not done this here since tendons do not contain many cells (cell fraction is measured at 2% by volume in adult rats [14] ), and to minimise the number of manual steps.Masking cells, however, would improve speed and efficiency as the FTA would not be seeking matches for non-fibrillar objects in the interfibrillar spaces.

Tracking accuracy
Fig. 3 B and D, together with video 1, illustrate the accuracy of the FTA in tracking fibres, in spite of the artefacts typically caused by this imaging modality.Small diameter fibrils are not captured as well as the larger diameter fibrils ( Fig. 3 B), the proportion of fibrils that are smaller than 150 nm that are tracked represents approximately 60% of the objects of this size seen in-plane.If this sample size is insufficient for a given study, it would be straightforward to trace the small numbers of missing smaller fibrils manually after running the FTA.This hybrid approach would still save the user time and resources, and enable them to capture a much greater field of view than is practical using manual segmentation alone.
To improve the automatic detection of small fibrils, however, examining short tracked segments which are not long enough for the minimum length threshold (see Section 2.5.1 ) to stitch them together could be a useful strategy.This would produce longer tracks and prevent them being discarded.To match short, jointed tracks, the error threshold could be lowered.Furthermore, small fibrils are more prone to being lost in the image segmentation step, so it could be beneficial to search more than one plane ahead in order to find the next segment.Alternatively, after running the whole FTA, the fibrils which are already successfully tracked could be stored and removed from the images, and the tracking algorithm could be run again using a different combination of α, β and γ in the error function.This would allow the capture of segments left behind in the first iteration of the algorithm.The latter option would be relatively easy to implement and could be repeated until no new tracks of the desired length are obtained.
Any of these techniques could be implemented in a research study where the small fibrils are of interest.For example, Szczesny et al. have suggested that small fibrils play a role in interfibrillar load transfer [73] , and Chang et al. [19] have found that the small fibrils are replenished throughout a circadian (24 h) time period.Though the capture of small fibrils could be improved in a number of ways, the FTA is a marked improvement on manual segmentation, which is commonly used [38,67] , and captures much of the tendon structure, as the larger diameter fibrils make up most of the volume.The FTA opens up a broader scope for microstructural and statistical analysis after tracking is complete, as much larger numbers of fibres can be traced than is practical manually.
Combination of the FTA with other computer vision techniques could be of benefit to improve tracking of small fibres.Almutairi et al. studied a similar problem using a mixed approach of machine learning and statistical techniques [74] .A supervised learning technique called Random Forests was used by Almutairi et al. to identify and classify fibres and a Kalman Filter was employed in the tracking stage to identify and later check and correct fibre connections.Classifying and correcting mistakenly disconnected short fibril tracks found by the FTA, using machine learning techniques, could improve its accuracy in future work.
Another potential improvement to the FTA could come from differential tracking, where the centroid error ξ c would not be calculated using the location of fibril i in plane p to compare to the centroid of fibrils in plane p + 1 ( Eq. ( 2) ), but a prediction of the location of the next segment in plane p + 1 would be made using a projection based on the recent history of the fibril's centroid movement between planes p − 1 and p.This technique comprises part of the work of Almutairi et al. discussed previously [74] .
We have tested implementing this approach in the FTA; however, when it encountered a damaged plane (discussed above, Fig. 3 C), the predicted centroid location was often erroneous, leading to the algorithm seeking the next segment in the wrong location.In this work, predictive mapping proved to be counter productive, as the algorithm was no longer flexible enough to overcome sectioning artefacts.This type of artefact is relatively common in SBF-SEM data, but would not be an issue for other 3D imaging techniques, such as XCT.Thus, it could be worth re-introducing this differential tracking method for data without sectioning artefacts, as it may be beneficial for fibres with weaker alignment with the longitudinal axis of the image stack than those in this study.
Of course, these additional potential features would add a computational cost.Presently, the algorithm takes around 2 h to run on a single 16 gigabyte core for an image stack of 700 image planes of 10 0 0 × 10 0 0 pixels.This is fast in research computing terms, and so this computational cost will not be a hindrance for future work.The image segmentation and tracking scripts do not currently make use of parallelisation, though for bigger jobs it may be worthwhile to do so to speed up computation.

Unusual structures: Branching, fibril ends, hairpins, rogue fibrils
Several unusual fibril features have been reported in the literature.Fibril ends are of interest in developmental studies [75] , and though rare across the depths typically covered by SBF-SEM, tapered fibril ends have been observed in tendon volumes [67] .It has been reported that fibril ends taper at a rate of approximately 30 molecules per 67 nm D-period [76] , and so the total tapering region would stretch over a few microns depending on the size of the fibril [21] .As the FTA does not require fibrils to have the same diameter along their entire length, but instead only requires approximate agreement locally (between adjacent planes, which are typically spaced 70-100 nm apart in SBF-SEM), the FTA would likely be able to track tapering fibrils.To identify the locations of fibril ends, a natural point to start is with the termini indicated in Fig. 3 D.In principle, these could either be true ends or fibrils that have not been successful matched to the next plane; however, since fewer than ten ends is typical in each plane, they can be inspected manually.We did not find any genuine fibril ends in the volume we analysed.
Highly tortuous small diameter fibrils often move much greater distances in-plane than their larger counterparts, often twisting and wrapping around larger diameter fibrils [19,67] .To account for this, the weighting of the centroid error in Eq. ( 2) could be made to vary with fibril diameter, allowing for greater lateral movement between planes for smaller fibrils.Hairpin-shaped fibrils (fibrils that turn back on themselves) have also been observed [21] .As the FTA typically traces fibrils from one end of the volume to the other, it cannot currently track loops of this kind.Hairpin bends will appear as two adjacent fibril tracks which merge and then disappear entirely.Similarly, fibril branching, where fibrils merge or separate, has been observed in tendons [21,40,67] .When fibrils branch, the algorithm will interpret this as one continuous fibril along one of the branches, and the other as terminating at the branch point.This is illustrated in Fig. 5 , where the orange fibril track splits into two tracks.Currently, the algorithm does not record branch points or hairpin bends and implicitly assumes fibrils are all continuous and are not connected; however, the presence of branched and hairpin-shaped fibrils in tendons could have an effect on tendon mechanical properties and fibril recruitment [21] .This is something that could be considered in the future, by introducing a criterion based on the centroids of the tracked fibrils becoming close to each other in regions in the vicinity of where one or both of the tracks end.

Conclusions
The FTA presented here can trace and characterise the geometry of large volumes of fibres from 3D imaging data, quantifying both the transverse and longitudinal geometry of the fibre population.We have used the algorithm to characterise the 3D microstructure of tendon fibrils, which are typically challenging to image, and their arrangement is difficult to characterise.This is due to their diameters being below the optical diffraction limit and their dense packing.In the field of tendon research, the FTA has the potential to be used in studies which utilise animal models and EM to study ageing, disease, development and injury.The FTA can enable this type of biological research to move beyond traditional 2D SEM to characterise and analyse microstructural features for different tendon phenotypes.We have developed the FTA in a general way, so it is suitable for analysing any material with approximately aligned fibres, such as blood vessels or ligaments, as well as engineered composite materials and textiles.In summary, the FTA is a powerful and flexible tool for studying structure-function relationships in a range of materials, is written exclusively using open source software, and is freely available to download and use.

Fig. 1 .
Fig. 1.Algorithm design and image segmentation.A -B : A randomly selected plane from a stack of SBF-SEM images presented in its original greyscale ( A ), and after the segmentation protocol outlined in Section 2.2 ( B ).The extent of the region of interest shown is 10.2 μm in both x and y .A small amount of cellular material can be seen in the top right corner.C : Schematic illustrating how the values of ξ (p , i, j ) ( Eq. (1) ) are generated for a pair of planes.For this simplified example of three fibrils in planes p and p + 1 , the 3 × 3 matrix of the values of ξ (p, i, j) would contain smaller values along the leading diagonal than the off-diagonal elements, as these are the correct (i, j) pairings.D : The red fibril segment i highlighted in plane p is compared to all the segments j in plane p + 1 in the white square, including those touching the border, and each value of ξ (p, i, j) is calculated and added to the error function matrix.This is repeated for every fibril segment i in plane p to evaluate the best matches for the

Fig. 2 .
Fig. 2. Results.A : 3D visualisation of tracked fibrils.B: Distribution of fibril lengths relative to the equivalent fascicle length.The red line is a log-normal probability density function fitted to the data.C: The minimum Feret diameter of the tracked fibril population.A multimodal distribution of fibril diameters is seen here, in agreement with previous work [19,68] .The red line is a bimodal normal probability density function fitted to the data.D : The orientation of the fibrils relative to the fascicle.The red line is a log-normal probability density function fitted to this data.E : The chirality of the fibrils, calculated using Eq.(13) .The dashed, black line is a skew normal probability

Fig. 3 .
Fig. 3.Model validation A : FTA accuracy at the plane in the middle of the volume.Red objects are tracked, grey objects are not tracked.B : The MFD of the tracked fibrils (blue) and the whole fibril population in a randomly selected 2D plane (grey).The dashed, black line represents the proportion of fibrils of each size that were successfully tracked, which ranges from approximately 60% to 100%.The FTA successfully tracks fibrils with an MFD of 150-300 nm, but is less successful in tracking fibrils with an MFD of less than 150 nm.C : A typical broken plane.This is seen where the sectioning is unsuccessful.In a volume of 700 planes, there are typically fewer than 10 broken planes like this.It is important that the FTA correctly handles these artefacts, otherwise erroneous matches or discontinuities in the fibril tracks will occur.D : The response of the FTA to broken planes.The arrows indicate the locations of broken planes within the volume, and the vertical axis indicates the number of fibril termini which are not close to the edge of the image (they are more than 5% of the width of the image from the border to exclude fibrils that leave the field of view from the calculation).There are typically fewer than 10 termini per plane, which, in principle, could either be true fibril ends or an artefact of tracking.Broken planes do not tend to lead to a large increase in termini, indicating that the FTA is handling them successfully.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Parameter optimisation in the error function.A : Illustration of the error threshold.The dashed line at ξ = 1 represents the error threshold used in this work.Each coloured line corresponds to a different, randomly chosen fibril, with a range of MFDs represented (see legend).The plot illustrates the 15 closest matches of each fibril segment in plane p to segments in plane p + 1 .All fibrils have a suitable match with the exception of that which has an MFD of 61 nm (represented by the darker blue line).B : Optimisation of the parameter values α, β, γ .The FTA was run with α = 1 , 0 ≤ β ≤ 3 and 0 ≤ γ ≤ 3 , and the proportion of fibrils that terminated in any given plane is plotted.A lower value (see legend) corresponds to a greater number of long fibril tracks, which are represented by darker blue areas in this heatmap.Partial fibril mappings for 50 random fibrils are plotted for three points of interest C : (i) represents the default values, (α, β, γ ) = (1 , 1 , 1) , (ii) and (iii) respectively correspond to when the value of α is high and low relative to β and γ .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Fibril branch point.The orange fibril branches to become two separate fibrils.The algorithm currently records them as two separate fibrils and does not record the branching.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Image processing pipeline for segmenting tendon fibrils from the extrafibrillar matrix and cells in SBF-SEM volumes.The original image, of size 4096 × 4096 pixels, was cropped to 1000 × 1000 pixels, which corresponds to 10.2 × 10.2 μm.