1.

Introduction

Second harmonic generation (SHG) imaging is a minimally invasive optical imaging technique with intrinsic optical sectioning capability.^1–3 The coherent nature of SHG leads to an intensity that is sensitive to the noncentrosymmetric arrangements of molecules⁴ and in biomedical applications, SHG has been shown to be effective in imaging fibrous collagen-rich tissues such as tendons, cartilage, and cornea.^5,6 Specifically, SHG microscopy has been successfully used in the label-free identification of structural changes in connective tissues under conditions such as skin photoaging, corneal abnormalities (keratoconus), liver fibrosis, and Osteogenesis imperfecta.⁷ In tissue engineering applications, SHG imaging has been useful in the quantification and monitoring of collagen growths in engineered cartilage tissues.^8–11 However, most efforts in using SHG imaging have focused on morphological characterization rather than quantification of the structural motifs of the fibrous structures. Therefore, combining SHG imaging with the appropriate methodology to quantify fiber structure can lead to objective evaluation in their structural organization. Recently, two-dimensional discrete Fourier transform (2D-FFT) has been shown to be an effective tool in analyzing fibril organization.^12,13 In particular, the average orientation and anisotropy of fiber tissue within any given region of interest (ROI) can be determined using 2D-FFT transformation. Both of these parameters vary with the domain size of analysis and the variation can be viewed as a quantitative standard to monitor changes in the fibrous tissues. However, in previous studies, the sizes of the imaging area used for 2D-FFT analysis are usually subjectively chosen.¹³ Therefore, the major purpose of the present study was to introduce a systematic approach to investigate the effect of the size of the analyzed region on the Fourier transformation analysis. In addition to local fiber orientation and aspect ratio, we introduce a quantitative parameter evaluating the correlation between the orientations of neighboring domains. For this purpose, the Ising-like order parameter was calculated for the fiber angular matrices. In addition to normal connective tissues, we show that our approach can be used to differentiate between normal and pathological fiber orientations in the cornea.

2.

Materials and Methods

2.4.

Fiber Orientation Analysis

To perform fiber organization analysis, the images were first divided into equally sized domains. Next, 2D-FFT was performed on each domain. To demonstrate our approach, we first applied our method on an artificially created fibrous pattern ($192\times 192\text{pixels}$) composed of adjacent domains, each $64\times 64\text{pixels}$ in size, with mutually orthogonal fibrous patterns in adjacent domains [Fig. 1(a) to 1(d)]. Next, the pattern was arbitrarily subdivided into domains of various sizes of $32\times 32$, $48\times 48$, $64\times 64$, and $96\times 96\text{pixels}$ in Fig. 1(a) to 1(d), respectively. The 2D-FFT analysis was applied on each subdivided domain, and the corresponding transformed images in the Fourier domain are shown in Fig. 1(a*) to 1(d*) in green. This procedure then was followed by a single pixel erosion-dilation process in order to smooth and delineate the contours of the transformed image.¹⁴ Quantification of the fiber geometry was then performed by fitting by an ellipse [shown in red in Fig. 1(a*) to 1(d*)] which allows the characterization of each domain by two parameters: (1) the angle between horizontal axis and the short axis of the ellipse, and (2) the aspect ratio (AR) between the short and long axes of the ellipse. The orientation of the short axis of the ellipse for each domain is shown by the white lines in Fig. 1(a) to 1(d). Since it is known that the short axis of the fitted ellipse is parallel to the average fiber orientation within each domain and that the AR describes the directional anisotropy, orientation properties associated with the fibers in each subdivided domain can be determined.¹³

Fig. 1

The use of different domain sizes for 2D-FFT analysis of an artificially produced pattern image. The images on the left show the division of the pattern image into domain sizes of (a) $32\times 32\text{pixels}$, (b) $48\times 48\text{pixels}$, (c) $64\times 64\text{pixels}$, and (d) $96\times 96\text{pixels}$. The corresponding Fourier transforms of the each region are, respectively shown in right along with the fitted ellipses [(a*) to (d*)]. The direction of the short axis of the fitted ellipse is shown by white lines for each domain in (a) to (d). The calculated average AR is also shown on corresponding images.

In this example, when fibers within each domain lie along the same direction, the 2D-FFT analysis accurately reflects the direction of the fiber orientation [Fig. 1(a) and 1(c)]. However, when the domain contains fibers along two orthogonal directions, the FFT analysis then reflects the average orientation of the fibers within the domain. Note, that in the latter case AR values are comparatively bigger [Fig. 1(b) and 1(d)]. According to our estimate, for a good correlation between the pattern line directions on images and those estimated from 2D-FFT the average AR value should be less than 0.5. Note that this criterion, as well as the optimal domain sizes, would depend on the characteristic extent of correlation in organization and needs to be derived for each case. In the case of collagen fiber factors, such as tissue type, image quality, resolution, and contrast would matter. Once these parameters are fixed, the corresponding estimate can be used unambiguously henceforth.

3.

Results and Discussion

Following a demonstration of our method with an arbitrarily generated pattern in Fig. 1, we applied our approach to more complicated fibrous tissues acquired through SHG imaging. Shown in Fig. 2(a) to 2(d) are the SHG images of RTT, MS, BC, and NHC, respectively. Next, the corresponding images were subdivided into different domains (each $32\times 32\text{pixels}$, $6.6\times 6.6{\mu \mathrm{m}}^2$) in which 2D-FFT was performed [Fig. 2(a*) to 2(d*)]. In these examples, since the fibers are not perfectly aligned along a given direction within each domain, the AR values derived from the fitted ellipses of the FFT image would increase compared to the case of perfectly aligned fibers. Also shown in Fig. 2(a*) to 2(d*) are red lines corresponding to the derived average fiber orientation within each domain.

Fig. 2

SHG intensity images of (a) RTT, (b) MS, (c) BC, and (d) NHC. Shown in (a*) to (d*) are the corresponding 2D-FFT analysis within the chosen domains (each $32\times 32\text{pixels}$ or $6.6\times 6.6{\mu \mathrm{m}}^2$). The bar on images is 40 μm.

As in the case of the arbitrarily generated fibrous pattern, the result of the 2D-FFT analysis is highly dependent on the size of the domain chosen. With a reasonably small domain size, the local fiber orientation can be better determined. However, a domain size too small would contain an insufficient amount of information of fiber morphology for accurate determination of fiber direction. With a large domain size, the 2D-FFT analysis would yield an average orientation within the domain and details of local fiber orientation would be lost.

To demonstrate these features, 2D-FFT analysis was performed using a range of domain sizes from $16\times 16$ to $512\times 512\text{pixels}$ which corresponded to actual image size of $3\times 3{\mu \mathrm{m}}^2$ to $100\times 100{\mu \mathrm{m}}^2$, respectively. Following the analysis, the average AR values were obtained from all domain sizes and the results are shown in Fig 3. In the cases of MS and RTT, the AR values are the result of averaging 10 images taken from 2 specimens; 5 from each. In the case of BC, ARs from images obtained at 10 different depths at 100 μm apart are averaged. NHC data represent values from forward SHG images obtained at 6 depths, 100 μm apart obtained from 2 specimens. A two-sided student’s t-test was used to compare parameters obtained for different specimens and the level of statistical significance was calculated using IBM SPSS software (Chicago, IL).

Fig. 3

Dependence of averaged ARs on the domain size chosen for 2D-FFT analysis of SHG images acquired from MS, RTT, BC, and NHC. Average ARs for MS and RTT are the results obtained from 5 images; averages from 10 images were calculated for BC and NHC. SHG images of corneas were obtained in the forward direction.

As Fig. 3 shows, the minima of the averaged AR values occur when the domain sizes are between $19\times 19$ and $36\times 36\text{pixels}$ (corresponding to $3.7\times 3.7$ and $8.2\times 8.2{\mu \mathrm{m}}^2$). Furthermore, the ARs of different fibrous architectures show similar trends at domain sizes below that at minimum AR. This observation is indicative of the fact that when domain size is decreased, the determination of general fiber orientation becomes erroneous. On the other hand, increasing the domain sizes will result in the inclusion of too many details of local fiber orientations, hence the averaging and loss of information on fine structure. However, as domain size increases beyond the point of minimum AR, our analysis starts to reveal differences between the different tissue morphologies. First, note that for the case of RTT, a relatively low AR value ($<0.4$) was obtained for domain size $96\times 96\text{pixels}$ ($19.8\times 19.8{\mu \mathrm{m}}^2$). This is indicative of the fact that collagen fibers in RTT specimens are well aligned. Note that at large domain sizes, the ARs for the four tissue types approached different terminal values of $0.40\pm 0.14$ (RTT), $0.76\pm 0.12$ (MS), $0.52\pm 0.08$ (BC), and $0.59\pm 0.03$ (NHC). Our results indicate the different collagen-containing tissues possess different degrees of short- and long-term order in organization.

Based on derived directional maps, we determined the Ising-like model correlation parameter according to Eq. (3). The average correlation parameter, calculated for different tissue types using the domain size at maximal AR value for the different tissues are presented in Table 1. Our results show that the calculated correlation parameters are consistent with the overall image morphology ranging from as low as $0.29\pm 0.08$ for the more random organization in MS dermis up to $0.75\pm 0.03$ for the highly organized RTT ($p<0.005$).

Table 1

Average correlation parameter of different tissue types.

As a final example of the potential applicability of our method in distinguishing normal and pathological tissues, we applied the 2D-FFT approach to analyze collagen fiber morphology in normal and pathological human corneas. Specifically, both forward SHG images from normal and keratoconical human corneas were obtained and analyzed. Representative large area images of normal and keratoconus cornea obtained at the depth of 100 μm are shown in Fig. 4(a) and 4(b). Slit-lamp examination of the keratoconus sample revealed a scar, which introduced additional turbidity into the orientation of collagen fibers of this sample [Fig. 4(b)]. The correlation parameter was calculated for entire images obtained at 6 depths for each specimen. AR values were calculated for different domain sizes and the minimum was determined to be $32\times 32$ ($6.6\times 6.6{\mu \mathrm{m}}^2$) for KHC and $36\times 36$ ($7.4\times 7.4{\mu \mathrm{m}}^2$) for NHC. Figure 4(a*) and 4(b*) show 2D-FFT analysis results of arbitrarily chosen ROIs from images 4(a) and 4(b) using the domain sizes corresponding to minimum average AR.

Fig. 4

Representative FSHG images of normal and keratoconus human corneas, as well as fiber orientation analysis of arbitrary selected ROIs from these images. (a) FSHG of NHC and (b) FSHG of KHC. (a*) and (*b) are results of 2D-FFT analysis of ROIs outlined in (a) and (b), respectively, analyzed using domain size of $36\times 36\text{pixels}$ for NHC and $32\times 32\text{pixels}$ for KHC. FSHG images were acquired at 100 μm depth.

In addition, the calculated average correlation parameters for the entire area of normal and keratoconus cornea images are shown in Fig. 5. The average correlation parameters corresponding to domain size $36\times 36\text{pixels}$ of minimum AR values are $0.53\pm 0.02$ and $0.43\pm 0.03$ for normal and keratoconus corneas, respectively ($p<0.005$). This serves as an example of the quantification of image morphology and shows that quantitative comparison of differences of fiber order can be performed.

Fig. 5

Average correlation parameter calculated using domain sizes of 26 to 48 pixel (5.3 to 9.8 μm) for NHC and KHC acquired by forward SHG at 100 μm depth. Arrows indicate the correlation parameter values corresponding to minimal AR values for both cases.

4.

Conclusion

In this study, 2D-fast Fourier transform analysis was developed and applied for morphological analysis of SHG images of different connective tissues. An approach, based on average AR estimation was developed, which allows for the evaluation of actual fiber directions. We found that both the aspect ratios and correlation parameter can be used to quantify differences in tissue morphology. While the aspect ratio analysis characterizes local fiber orientation within the ROI domains, the correlation parameter is effective in revealing the orientation order among the ROIs. To avoid ambiguity, our analysis was performed on systematically selected ROI sizes and we found that the dependence of AR and correlation parameters on ROI sizes can be used as additional metrics for characterization fiber orientation. In addition, to distinguish fiber morphologies in normal tissues such as RTT, MS, or BC, we found that organizational differences between normal and pathological cornea (keratoconus) can be achieved. The demonstrated approach can be used for automated morphological analysis of tissue fibers and may be extended for diagnostic purposes.