1.

Introduction

Defining fiber orientation can be a critical step in assessing the characteristics of a variety of engineered and native tissues. For example, assessing electrospun biomaterial fiber orientations is key for developing a composite biomaterial with desired anisotropic properties.¹ In fact, a variety of microstructural models of tissue biomechanics require fiber orientation measurements to accurately predict mechanical properties.^2,3 Additionally, fiber orientation distributions can be an important metric in diagnosing diseased tissues, localizing injuries, or detecting engineered tissue development.^4–8 To create maps of collagen fiber orientations within a tissue, techniques such as quantitative polarized light imaging, small angle light scattering, and polarization-sensitive second harmonic generation imaging have been employed.^9–12 However, these techniques require additional optical elements and/or hardware that are not frequently available in the laboratory. For the majority of research projects, images are acquired in which fiber orientation is not an inherent source of contrast. Rather, fiber orientation can be qualitatively identified by acquiring an image of sufficient magnification to spatially resolve individual fibers organized within a tissue. Depending on fiber size, imaging approaches can range from using a simple light microscope to second harmonic generation imaging to scanning electron microscopy.

A number of current approaches exist to quantify the average fiber direction in an image. Most commonly, a Fourier transform is applied and a two-dimensional (2-D) power spectral density (PSD) image can be computed.^13–16 Sampling the average power of the PSD at pixels corresponding to different angles can yield an accurate measure of the total fiber orientation distribution in the image.⁶ However, such techniques generally require a square image size and can be prone to systematic measurement errors produced by large scale features such as the absence of fibers in irregularly shaped subregions within the image. The Hough transform can also be utilized to identify fiber orientation information by characterizing the polar coordinates of specific image locations. By making cumulative measurements of the polar coordinates associated with weighted pixel locations, the Hough transform has been used to determine collagen fiber orientations within $32\times 32\text{pixel}$ subregions in SHG images.⁶ This image processing technique has also been applied to identify muscle fiber orientation from ultrasound images, and is amenable to any modality with sufficient contrast and resolution to visualize discrete fiber bundles.¹⁷ However, to generate pixel-wise fiber orientation maps, performing Fourier or Hough transforms surrounding each pixel within an image is particularly time consuming.

More sophisticated fiber tracking techniques have also been developed, in which fiber objects are defined and tracked through energy minimization (e.g., active contours) or line propagation algorithms.^18–20 However, these advanced techniques utilize iterative processes, which are relatively time consuming and typically require some level of user interaction. As a result of these limitations, fiber orientation data are not traditionally calculated or reported in microscopy studies. Yet, the ability to automatically generate fiber orientation information at each pixel within an image would be advantageous for disease diagnosis, assessing tissue development, and many mechanobiology applications.

A few simple computational approaches have been developed for rapidly detecting fiber orientation in micrographs. Traditional edge detection algorithms use Sobel or Canny operators to compute image gradients in the $x$- and $y$-directions, find local maximums in the magnitude of the gradient, and identify the location of fibers within an image.^21–23 In addition to identifying edges, these algorithms can estimate orientation by calculating the arctangent of the $x$ and $y$ gradient at each pixel. However, this approach generally suffers from inaccuracy that makes it unsuitable for quantifying exact fiber orientations. The goal of this study is to develop an algorithm for the rapid, yet accurate quantification of fiber orientation at each pixel within an image. A weighted orientation vector summation approach has been developed in MATLAB that can provide pixel-specific orientations at an order of magnitude faster than pixel-wise Fourier analysis with improved accuracy over gradient-based edge detection approaches. The accuracy of this detection algorithm is assessed by simulating images with known orientations, and computational time is monitored for a variety of image and kernel sizes. Applicability of this technique to a variety of imaging modalities is also demonstrated.

2.

Methods

2.1.

Overview of Vector Summation Technique for Orientation Detection

The algorithm developed in this study defines fiber orientation by identifying the variability of image intensities along different directions surrounding each pixel within an image. Upon selecting a square window size to evaluate each pixel, all possible vectors passing through the center pixel in the window and two additional pixel locations symmetric about the center are defined, and the angles associated with these orientation vectors are calculated [Fig. 1(a)]. Next, for each pixel location, the orientation vectors are weighted based on intensity variability and length [Fig. 1(b)]. Because the perimeter of a circle, with diameter $L$, is $\pi *L$, the lines of length $L$ passing through the center pixel of the window also scale in number by $\pi *L$. Since the longer vectors are more numerous, they would have a greater collective weight when calculating mean fiber orientation surrounding the center pixel. To correct for this effect, each orientation vector is weighted by the inverse of its length ($L$), which enables pixels close to the center pixel of interest to have equal collective weight in determining fiber orientation as pixels farther from the center [$W_1$ in Fig. 1(b)]. Sensitivity to fiber orientation was determined by also weighting the vectors according to the lack of variance among the intensities of the three pixel locations used to define the vector. This vector weight was defined as the difference of the maximum possible standard deviation among three points and the measured standard deviation among the points [$W_2$ in Fig. 1(b)]. The images were converted to span intensities from 0 to 1, and the maximum theoretical standard deviation among three values spanning such a range is the square root of $1/3$. Once the weighted vector lengths are computed for a given pixel, all the vectors are summed to determine an average direction [Fig. 1(c)]. While vector summation can provide mean orientation data of circular data (for which an orientation of 0 deg and 360 deg are equivalent), fiber orientations are axial data (for which 0 deg and 180 deg are equivalent). To account for these axial data, all vector angles are multiplied by two prior to vector summation and the angle of the resultant vector is divided by two in order to determine the final fiber orientation angle at the center pixel.

Fig. 1

Summary of vector summation technique for pixel-wise fiber orientation measurements. The vector lengths in (b) and (c) represent the relative weighting factors for each orientation. The average vector orientation corresponds to the direction with the lowest variation in pixel intensity values ($a_i$).

The time required to serially compute the fiber orientation for every pixel within an image using this vector summation method, or any other algorithm, would normally be prohibitive. However, this vector summation process can be performed simultaneously at each pixel to dramatically increase the speed by which an orientation map can be generated. The initial vector orientations and lengths can be computed once at the beginning of the algorithm. Next, our program creates three copies of the intensity image and, based on the Cartesian coordinates of a given orientation vector, shifts the coordinates of two copies in opposite directions. These three copies of the image with different coordinate shifts are assembled into an $n\times m\times 3$ array of intensities (where n and m are the image dimensions). For example, to analyze a $512\times 512$ image using a $7\times 7$ window there are 24 unique vectors [Fig. 1(a)], and for each one of those vectors, an array with dimensions of $512\times 512\times 3$ is produced. For a vector with coordinates of (2,1) relative to the center pixel (26.6 deg above the horizontal), one $512\times 512$ copy is shifted to the right 2 pixels and up 1 pixel, while another copy is shifted to the left 2 pixels and down 1 pixel. The standard deviation along the third dimension of this three-dimensional (3-D) array is computed to obtain the $W_2$ weight for a given orientation vector at each pixel. This process for determining the each vector’s weight is repeated for each vector position. Every time a weighted vector is computed, the $x$- and $y$- components of the vector are each added to running subtotals of the components for each pixel. The cumulative $x$- and $y$- component vectors are ultimately used to compute the average vector orientation at the center pixel. With this algorithm, the computational time for orientation detection more closely scales with the number of orientation vectors (i.e., the window size) used rather than the number of pixels in the total image.

3.

Results

The accuracy of the average fiber orientation measurement using the vector summation method depends on both the thickness of the fibers and the window size used in the analysis. When the proximity of fibers to each other is not a factor, a larger window size used for detection will produce smaller errors in orientation measurements regardless of fiber thickness (Fig. 2). For most window sizes, error in the fiber orientation measurements increases exponentially as the fiber thickness approaches, and eclipses the width of the window size (Fig. 2). Interestingly, for window sizes of $9\times 9$ and smaller, the detection algorithm is most accurate when the thickness of the simulated fibers is 5 pixels (Fig. 2). Although these results suggest larger window sizes will result in improved accuracy, these simulated images maintained inter-fiber spacing that was greater than the maximum window size evaluated.

Fig. 2

Error in fiber orientation measurements with respect to fiber thickness and window size used for the vector summation measurements. Beyond a thickness of 5 pixels, error in fiber orientation measurements began increasing exponentially when fiber thickness increased beyond the window size used for detection. Standard deviation of the average difference in angles was on the same order of magnitude as the mean value, but is not plotted in the graph in order to aid in viewing the overall effects of thickness and window size.

The accuracy of fiber orientation detection as a function of window size differs from Fig. 2 when the simulated images contain fibers of different orientations, thicknesses, and proximity to each other. Specifically, measurable improvements in the average accuracy are only detected between $3\times 3$ to $11\times 11$ window sizes (Fig. 3). Using an $11\times 11\text{pixel}$ window size, the absolute value of the difference in detected orientation is $2.5\pm 2.7\mathrm{deg}$. There is no measurable source of systematic error in the measurements as the average difference in angle is $0.00\pm 3.68\mathrm{deg}$. For window sizes greater than $11\times 11\text{pixels}$, negligible changes are identified in the average absolute error and a minimum mean absolute error of $2.2\pm 2.9\mathrm{deg}$ is reached using a $15\times 15$ window. As the window size increases, the standard deviation of the average error values increase from a minimum of 2.7 deg using a $13\times 13$ window to a maximum of 5.7 deg for a $25\times 25\text{pixel}$ size. While the error for the majority of pixels improves with an increase in window size, an increase in error is produced for certain locations if the window size is large enough to include multiple simulated fibers with different orientations in close proximity to each other [Fig. 3(d)]. This increased spatial variability in error values produces an increase in the standard deviation of angle differences without an increase in the mean difference [Fig. 3(a)].

Fig. 3

Error in orientation detection using the vector summation method for a simulated image containing regions with diverging orientations. (a) The average absolute value of the difference between measured and actual pixel-wise orientations noticeably decreases with increasing window sizes of $3\times 3$ to $11\times 11$ using a simulated image (b) of circles at different locations, diameters, and thicknesses. (c) With increasing window size, error initially decreases and then becomes localized to regions where fibers are in close proximity to each other. (d) A magnified view of the error in the area indicated by red outline in (b) demonstrates an increase in error between $11\times 11$ and $25\times 25\text{pixel}$ windows for specific regions with multiple nearby fibers.

When compared to a pixel-wise Fourier-based detection of fiber orientation, small differences in accuracy are identified, but computational time is an order of magnitude faster in the vector summation method. Specifically, a Fourier-based approach produces an average absolute difference in fiber orientation of $2.7\mathrm{deg}\pm 1.8\mathrm{deg}$ using an $11\times 11$ window, compared with $2.5\mathrm{deg}\pm 2.7\mathrm{deg}$ for the vector summation approach. Interestingly the error of the Fourier-based measurements could be further reduced to $1.4\mathrm{deg}\pm 1.2\mathrm{deg}$ by multiplying each pixel in the 2-D PSD by the distance from the center pixel, which effectively gives more weight to the smaller frequencies when calculating the mean fiber direction. Although the vector summation technique performs with slightly lower accuracy than a frequency-weighted Fourier approach, the time required to detect alignment at each pixel in a $1920\times 1920$ image using an $11\times 11$ window is 10-fold shorter. Computational time scales exponentially with image size for both detection methods [Fig. 4(a)], but the vector summation method becomes increasingly advantageous with larger image sizes. Similar relative increases in computational time are observed between methods for increasing window sizes [Fig. 4(b)].

Fig. 4

Comparison of the computational times for a pixel-wise Fourier detection of fiber orientation and the vector summation approach. (a) The vector summation approach ranges from 5.25-fold faster for a 0.01 megapixel image to 10.34-fold faster for a 5 megapixel image compared to Fourier approaches for a range of image sizes using an $11\times 11$ window size. (b) Smaller window sizes produce a more pronounced difference in the computational time of the two methods. Using a $960\times 960$ image, computational time of the vector summation method was 74.4-fold faster using a $3\times 3$ window, but reduced to a 3.2-fold difference using a $25\times 25$ window.

Sobel and Canny-based directional gradient detection algorithms provide pixel-wise orientation maps 1 to 2 orders of magnitude faster than the vector summation approach developed in this paper, but the accuracy of the directional gradient algorithms is also far worse. Using a $3\times 3$ Sobel filter, the mean absolute difference in orientation measurements of the image in Fig. 3(b) is $12.1\mathrm{deg}\pm 18.8\mathrm{deg}$. A more complex $10\times 10$ Canny operator provides a mean absolute difference of $9.9\mathrm{deg}\pm 16.6\mathrm{deg}$. No measurable improvements in accuracy were detected by varying the kernel size of the Canny operator, demonstrating the vector summation technique developed here can provide fiber orientation detection with a 4-fold increase in accuracy.

To assess the broad applicability of this approach, fiber alignment distributions were calculated from biomedical images that were acquired using a range of modalities (Fig. 5). The isotropic alignment of neurites sprouting from a chick dorsal root ganglion (DRG) explant body was measured from a standard fluorescence microscope [Fig. 5(a)]. Electrospun silk fibers were imaged using scanning electron microscopy (SEM) and the strong alignment towards 135 deg was evident from the orientation maps and distribution histogram [Fig. 5(b)]. Although a SHG image of collagen fiber deposition in engineered bone tissue demonstrated less anisotropy than electrospun silk fibers, a greater proportion of fibers near 45 deg was evident from the false colored orientation map and the pixel-wise distribution histogram [Fig. 5(c)]. Light transmission images of MCF10A cells within a collagen matrix demonstrated a preferred orientation towards 120 deg [Fig. 5(d)]. With each image modality, the calculated pixel-wise orientation measurements matched the qualitative observations of fiber alignment from the original image.

Fig. 5

Fiber orientation detection via vector summation performed on (a) fluorescence microscopy image of DRG neurites, (b) scanning electron microscopy image of silk fibers, (c) second harmonic generation image of collagen in engineered bone, (d) picosirius red stained collagen fibers surrounding MCF10A cells (reproduced with permission²⁴). For the orientation maps in the middle row, pixel-wise fiber angle measurements were color coded to the HSV color map in MATLAB, and these color maps were then multiplied by the grayscale intensity image to produce a final image.

4.

Discussion

A weighted vector summation method for detecting pixel-wise fiber orientations from medical images has been developed. This technique relies on changes in the fiber intensity variations in different directions and can provide orientation information with 2-3 deg mean accuracy an order of magnitude faster than the evaluation of orientation using a pixel-wise Fourier-based approach (Fig. 4). The enhanced speed of this algorithm relative to Fourier approaches now makes the definition of pixel-wise fiber orientation maps practical for experiments that generate large series of images or particularly high resolution images. Although pixel-wise fiber orientations can be estimated at an additional 1 to 2 orders of magnitude faster by computing directional gradients with a Canny operator^21,23 (approximately 1.2 s in this study), the large error in these measurements makes such a technique insufficient for most quantitative analyses of biomedical images. The improved accuracy of the developed vector summation technique over Canny-based approaches and the improved computational speed over Fourier-based pixel-wise measurements enable this technique to be uniquely amenable to quantifying and visualizing different fiber orientations for a range of imaging modalities and applications.

When implementing this technique for detecting fiber orientation, a window size that is greater than the diameter of the fibers should be selected (Fig. 2). Interestingly, our analysis of simulated fiber images indicates that the detection algorithm is most accurate when the fiber diameter is at least 5 pixels, suggesting the digitization of the image may reduce the accuracy of orientation measurements of fine fibers (Fig. 2). Selecting a much larger window size than the average fiber diameter will not necessarily produce more accurate measurements when fibers of different orientations are in close proximity to each other (Fig. 3). This is because a larger window size essentially blurs the local orientation with larger surrounding areas. Additionally, when the inherent contrast of the raw image produces intensity variations across the diameter of the fiber, such as in SEM images [Fig. 5(b)], this method appears particularly sensitive, regardless of fiber diameter. Furthermore, because fiber orientations are governed by directional sensitivity to intensity variance and not mean intensity, accurate orientation information is not limited to locations along the center-line of the fiber [Fig. 5(b)].

The ability to quantify fiber orientation at each pixel through this technique comes at the cost of some precision by which an accurate orientation distribution can be determined. Traditional Fourier analysis transforms the entire image to Fourier space, providing power densities for a much greater range of frequencies and orientations than a small $11\times 11$ subimage. As a result, the fiber orientation distribution can be defined with greater precision than any pixel-wise detection technique. For applications that require the definition of fiber orientation with subdegree precision, with no need for an understanding of the spatial distribution of the orientations within the image, calculating the power spectral density via a Fourier transform is still most appropriate. However, quite often medical images contain additional features such as cells or tissue edges that can cause systematic errors in Fourier-based approaches. For images with nonfiber features, such as the elongated cell in Fig. 5(b), a bias in the measured orientation would be produced due to the alignment of the cell. Although the vector summation technique proposed in this study also produces a 45 deg artifact within the cell region in Fig. 5(d), a simple mask to remove these pixels from the analysis can be implemented. Masking the cell regions prior to a Fourier transform, however, will not necessarily remove orientation measurement bias, because the absence of signal in the shape of a cell would also produce a bias in the direction of the cell orientation. In summary, while lacking the subdegree precision of approaches using whole image Fourier or Hough transforms, pixel-wise techniques can more easily account for systematic errors produced by large scale image shapes and features by masking image regions in the spatial domain.

In the algorithm presented in this study, a set of vector orientations [Fig. 1(a)] within a specific window size are weighted and summed to compute an average orientation at each pixel. Rather than serially computing the orientation at each pixel, running totals of the $x$- and $y$- components of the weighted orientation vectors are computed for all pixels simultaneously within an image in MATLAB. By serially summing the different vectors within the window in all pixels simultaneously, rather than serially performing the vector summation at each pixel, computational time does not scale up with image size as rapidly as a pixel-wise Fourier approach [Fig. 4(a)]. Although Fourier-based approaches can utilize an analogous vector summation approach to compute the mean orientation, a discrete Fourier transform must be applied surrounding each pixel prior to vector summation, which results in the substantial increase in computational time (Fig. 4). The ability to remain in the spatial domain provides a clear advantage in the computational time for this vector summation technique, but parallel processing techniques may be used in the future to overcome some of the limitations in Fourier approaches.

The pixel-wise weighted vector summation technique developed in this study provides superior accuracy ($2.5\mathrm{deg}\pm 2.7\mathrm{deg}$ mean error) in orientation measurements compared to pixel-wise measurements from gradient-based filters using a standard Canny operator ($9.9\mathrm{deg}\pm 16.6\mathrm{deg}$ mean error). Although a small improvement in accuracy can be accomplished by using a pixel-wise frequency-weighted PSD technique ($1.4\mathrm{deg}\pm 1.2\mathrm{deg}$ mean error), the weighted vector summation technique provides orientation maps an order of magnitude faster than any Fourier approach. As a result, the method developed here makes the real-time production of accurate pixel-wise fiber orientation maps possible for machine vision applications, real-time medical diagnostic tools, and basic science research. The ability to produce pixel-wise orientation data can enable a wider variety of diagnostic metrics beyond directional statistics based on the total orientation distribution of the image that is normally produced through traditional analyses. Specifically, directional summary statistics can be computed from irregular shaped areas, the spatial autocorrelation of directional variation can be determined to assess fiber organization, and unique measures of fiber undulation (i.e., crimping) can be produced. Additionally, the definition of pixel-wise fiber orientation maps can help facilitate the implementation of advanced fiber tracking techniques that rely on directional information provided by modalities such as diffusion tensor MRI or polarization microscopy.¹⁸ In summary, the speed by which fiber orientation maps can be generated using this vector summation technique enables the practical implementation of unique directional outcomes for most experimental and clinical applications in which tissue microstructure is imaged.