## 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.^{1} 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}5.6.7.^{–}^{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}10.11.^{–}^{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}14.15.^{–}^{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.^{6} 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{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixel}$ subregions in SHG images.^{6} 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.^{17} 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}19.^{–}^{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}22.^{–}^{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.

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.

## 2.2.

### Evaluation of the Vector Summation Detection Technique Performance

To evaluate the accuracy of this vector summation method to detect fiber direction, images with known orientations were created and evaluated. In an effort to include all possible fiber orientations in each image, white circular rings were created over a black background. First, to determine the relationship between detection window size and fiber thickness, three rings of 100, 200, and 300 pixel radii were produced. Binary images of these rings were produced with thicknesses of 1, 2, 3, 5, 7, 11, 15, and 25 pixels. A $9\times 9$ Gaussian filter with $\sigma =1.2$ was applied to the images to smooth the rings and create a grayscale image. The orientation at each pixel contained within the ring was defined by automatically detecting the center of each ring and finding an angle that is orthogonal to a vector from the center of the circle to the pixel of interest. To further investigate the accuracy of this technique, an additional image with circular rings of varying diameter, thickness, and orientation was created with circles positioned at varying distances from each other. The absolute value of the difference between the true orientation and the measured orientation was calculated for each pixel, and this average absolute difference in angles was computed for measurements using a variety of window sizes ranging from $3\times 3$ to $25\times 25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$.

To provide context for the error in the vector summation-based measurements, the absolute difference in angle measurements using this technique was compared with a pixel-wise Fourier-based approach. Based on previously described methods, fiber orientation was calculated from sampling different angles in a 2-D power spectral density (PSD) image.^{6} Unlike previous work, this PSD calculation was performed at each pixel on subimages of the same window size as the vector summation technique. Prior to transforming each subimage to Fourier space, the intensities were apodized by multiplying by ($1/r-1/{r}_{\mathrm{max}}$), where $r$ is the distance from the center pixel and ${r}_{\mathrm{max}}$ is the radius of the largest circle that could fit into the subimage space. This apodization darkened the edge of the subimage and insured a bias toward 0 deg or 90 deg would not be present due to the discontinuities between the opposite edges of the subimage. Traditionally, once the 2-D PSD is computed, the power density is averaged over discrete angular segments to determine the orientation distribution and ultimately the mean fiber direction.^{6}^{,}^{14} Due to the limited data within the subimages of sizes $25\times 25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$ and smaller in this study, a summation of PSD-weighted orientation vectors corresponding to each pixel location was performed to calculate the mean fiber direction. The accuracy and computational time for a 2.53 GHz Intel Core i5-460M processor using a range of window sizes and square image sizes was determined and compared to the spatial domain vector summation technique developed in this study.

In addition to comparisons with a pixel-wise Fourier technique, pixel-wise orientations were calculated using directional gradients derived from Sobel or Canny operators.^{22}^{,}^{23} To generate an image gradient in the $x$-direction using a Sobel operator, a $3\times 3$ kernel was produced by multiplying ${(\begin{array}{ccc}1& 2& 1\end{array})}^{T}$ by $(\begin{array}{ccc}-1& 0& 1\end{array})$ and then was convolved with the image. This process was repeated for the y-direction by taking the transpose of the $3\times 3$ kernel and performing another convolution. Fiber direction was inferred by calculating the arctangent of the $x$- and $y$- gradients using the atan2 function. In addition to a Sobel-based detection, this study investigated the accuracy of fiber directions derived from Canny operators as well. A $10\times 10$ Canny operator was defined by computing the first derivative of a Gaussian distribution in either the $x$- or $y$-direction, and convolved with the image. The Canny and Sobel algorithms were applied to the simulated image containing circles with variable properties and positions, and as with the other techniques, error was defined by computing the average absolute difference in the actual and measured fiber orientations.

## 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.

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{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixel}$ window size, the absolute value of the difference in detected orientation is $2.5\pm 2.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. For window sizes greater than $11\times 11\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$, negligible changes are identified in the average absolute error and a minimum mean absolute error of $2.2\pm 2.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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{\hspace{0.17em}}\text{\hspace{0.17em}}\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)].

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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 1.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ using an $11\times 11$ window, compared with $2.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 2.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ for the vector summation approach. Interestingly the error of the Fourier-based measurements could be further reduced to $1.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 1.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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)].

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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 18.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. A more complex $10\times 10$ Canny operator provides a mean absolute difference of $9.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 16.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.

## 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm \phantom{\rule{0ex}{0ex}}2.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ mean error) in orientation measurements compared to pixel-wise measurements from gradient-based filters using a standard Canny operator ($9.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 16.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ mean error). Although a small improvement in accuracy can be accomplished by using a pixel-wise frequency-weighted PSD technique ($1.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}\pm 1.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.^{18} 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.

## Acknowledgments

This research was supported by Grant Numbers R01EB007542 from NIBIB/NIH and F32AR061933 from NIAMS/NIH. We would like to thank Amy M. Hopkins [Fig. 5(a)], Rodrigo R. Jose [Fig. 5(b)], Dr. David L. Kaplan [Fig. 5(a) and 5(b)], Dr. Carlo Amadeo Alonzo [Fig. 5(c)], Dr. Ana M. Soto [Fig. 5(d)], Dr. Carlos Sonnenschein [Fig. 5(d)], and Dr. Eugen Dhimolea [Fig. 5(d)] for providing images to demonstrate the capabilities of our technique.

## References

*In vivo*, pixel-resolution mapping of thick filaments’ orientation in nonfibrilar muscle using polarization-sensitive second harmonic generation microscopy,” J. Biomed. Opt. 14(1), 014001 (2009).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3059627 Google Scholar