Open Access
Add to BibSonomyAbstract
Optical polarization tractography (OPT) has recently been applied to map fiber organization in the heart, skeletal muscle, and arterial vessel wall with high resolution. The fiber orientation measured in OPT represents the 2D projected fiber angle in a plane that is perpendicular to the incident light. We report here a dual-angle extension of the OPT technology to measure the actual 3D fiber orientation in tissue. This method was first verified by imaging the murine extensor digitorum muscle placed at various known orientations in space. The accuracy of the method was further studied by analyzing the 3D fiber orientation of the mouse tibialis anterior muscle. Finally we showed that dual-angle OPT successfully revealed the unique 3D “arcade” fiber structure in the bovine articular cartilage.
© 2016 Optical Society of America
1. Introduction
Fibrous tissue exists in many parts of the body, for example in skeletal muscle, heart, and brain. A well-organized fiber architecture is essential for maintaining normal physiological functions in these tissues. As an example, the unique helical myofiber organization in the heart is critical for regulation of blood pumping [1]. A technology that can image fibrous organization with sufficient resolution and field-of-view would be valuable for our understanding and assessment of the important structure-function relationship in these tissues [2, 3]. Traditional histology sectioning [4] can achieve cellular level resolution, but is labor-intensive and only practical for sampling a small area in fixed tissues. Diffusion-tensor based magnetic resonance imaging (DTI) [5] has been established for imaging global 3D fiber organization, especially in brain and heart. However the spatial resolution of current DTI is limited to sub-millimeters which is insufficient for imaging thin tissues or small structural changes.
High resolution optical coherence tomography (OCT) has been investigated to study fiber orientation in tissues [6–10]. However, fiber detection based on intensity contrast can be affected by non-fiber related intensity variations and only works in tissues with large fibers. We recently developed optical polarization tractography (OPT) for imaging tissue fiber orientation [11–13]. OPT is based on Jones matrix polarization-sensitive OCT (PSOCT) and uses the depth-resolved “local” optic axis to construct the 3D fiber tractographic images of the sample. OPT has been recently applied to imaging fiber structure in heart muscle [14], skeletal muscle [15], and blood vessels [16]. Histology validation showed that OPT can image detailed fiber architectures in tissue with histology-like resolution [17].
The fiber orientation measured in OPT represents the “projected” orientation in an imaging plane perpendicular to the incident light [17]. This projection angle in a 2D plane cannot fully represent the true fiber orientation in the 3D space, especially when a fiber is inclined toward the incident light. We show in this study that images of true 3D fiber orientation can be constructed by applying OPT to image a sample at two different angles. Using more than one incident angle to resolve the absolute fiber orientation in polarization imaging has been explored before in a few studies [18–23]. For example, a “dual-projection” polarimetry imaging system was developed to obtain orientation map in tissue slice [18]. Analyzing cumulative PSOCT results obtained at different incident angles was also used to derive the representative 3D fiber orientation in cartilage [20–22] and muscle [23]. Unfortunately, these earlier studies did not reconstruct the depth-resolved 3D fiber orientation in the sample. In sharp contrast, here we show that dual-angle OPT successfully reveals the full high-resolution 3D images of 3D fiber orientation in complicated biological tissues such as the articular cartilage.
2. Methods
2.1 OPT algorithm and imaging system
The OPT method uses Jones calculus based algorithms to calculate the “local” optic axis from cumulative measurements in Jones-matrix PSOCT as described in detail previously [13,14]. Briefly, the “round-trip” pixel-wise Jones matrix of the entire imaging volume JRT(ϕn,2ρn) was constructed, where the cumulative complex retardance ρn and axis ϕn from the nth depth layer can be calculated. With ρn and ϕn, a new Jones matrix JL(ϕn, ρn) was constructed and another complex matrix K(n + 1) was generated as:
The “local” complex retardance at the nth layer ξL(n) = δL + iσL was obtained using the eigenvalues of K(n + 1). The imaginary σL was eliminated from JRT(n) in subsequent calculations to create a diattenuation-free matrix J′RT(n). The “local” optic axis β was derived using an iterative process from the first layer where the axis was obtained directly from the round-trip Jones matrix J′RT(1) using eigen-decomposition. The derived “local” axis and retardance from the 1st to the (n-1)th layer were used to construct the “single-trip” Jones matrix JST(n-1). Then the nth layer local axis was obtained from the local Jones matrix JL(n):The local optic axis (β) was used to construct the tractography using the streamline functions in the Matlab software.The OPT algorithms are general and can be applied to any Jones matrix PSOCT systems. For the purpose of demonstration, we used a single-camera spectral-domain Jones-matrix PSOCT in this study. The details of this system has been reported elsewhere [24,25]. It utilized a 0.85 μm superluminescent diode (SLD) as the light source and had an axial resolution of 8.2 μm in air and a lateral resolution of 12.4 μm. The system acquired a 6 × 6 × 1.1 mm3 sample volume (2000 × 1000 × 280 pixels in B × C × A scans) at a 50 kHz A-scan rate. In this study, the obtained 3D data set was interpolated in the lateral and axial directions to achieve an isotropic pixel size of 6.0 µm in all A-, B-, and C-scan directions.
2.2 Dual-angle OPT
Figure 1 shows the laboratory coordinate system established using the orthogonal A-, B-, and C-scan directions labeled as the z-, y-, and x-axis, respectively. The origin is located at the intersection of the incident light and the sample surface (Fig. 1). The incident light is refracted at the tissue surface depending on the incident angle and enters the tissue as the z’-axis. The local optic axis measured in OPT is equivalent to the orientation of the actual fiber projected at the “evaluation plane” (the x′-y′ plane in Fig. 1) which is determined based on the law of optical refraction in 3D. This evaluation plane is perpendicular to the refracted light (z′-axis) inside the tissue, the incident plane (z-z′-x′), and the projection plane (formed by z′ and the 3D fiber). It is important not to confuse the “evaluation plane” with the conventional en face plane (x-y plane) in OCT imaging. In addition, OPT measures fiber orientation (β) with respect to the horizontal linear polarization in our system, which is aligned with the x-axis in the laboratory coordinates. This reference direction also needs to be projected to the evaluation plane by rotating θ around the y′-axis (normal to the incident plane) and is marked as the “angle reference” (0°) in Fig. 1.
Fig. 1 An illustration of the OPT measurement geometry. The laboratory coordinate system is formed by the x-y-z vectors. The light is initially incident along the z-axis, but is refracted to the z′-axis inside the sample. OPT measures the fiber orientation β within the evaluation plane (x′-y′) which is perpendicular to the light direction (z′) inside the sample. The actual fiber lies within the projection plane formed by the light direction and the fiber.
The refracted light (the z′ axis in Fig. 1) can be obtained based on the Snell’s law using the incident light vector I = [Ix, Iy, Iz]T and the surface normal vector N = [Nx, Ny, Nz]T:
where n is the tissue refractive index. The θi is the incident angle so that cosθi = −N•I. The refracted light vector can also be represented as L = [sinθcosϕ, sinθsinϕ, cosθ]T where the polar angle θ and azimuthal angle ϕ are defined as θ = cos−1(|Lz|), ϕ = tan−1(Ly/ Lx).The orientation β measured in OPT defines a 2D unit vector which represents the projected fiber inside the evaluation plane (x′-y′ in Fig. 1). This projected vector (Fig. 1) can be expressed in the laboratory coordinate (x-y-z) as:
The actual 3D fiber cannot be determined by the β measured in a single-scan OPT measurement. However this fiber must be located within the projection plane (Fig. 1) formed by the refracted light and the projected fiber. If the fiber orientation is measured using a different projection plane, the 3D fiber orientation can be fully determined as the intersection line of the two projection planes. The second measurement can be performed by rotating the sample around a predefined axis. The projection plane can be conveniently identified using its normal vector M by the cross-product of L in Eq. (3) and P in Eq. (4):In the second measurement, the refracted light has a different direction (θ′, ϕ′) leading to a new projected fiber orientation β′ and projection plane with a normal vector M′(θ′, ϕ′, β′). The 3D fiber orientation vector F of the sample was calculate as:To facilitate the evaluation of 3D orientation, the constructed 3D fiber vector F at each pixel can be projected into the three orthogonal planes y-x (B-C), z-x (A-C) and z-y (A-B) to obtain the corresponding 2D projection angles:where θf and ϕf are the polar and azimuthal angles of the 3D fiber direction.2.3 Image registration
In this study, dual-angle OPT was implemented by imaging the sample placed at two different rotation angles (γ1 and γ2) around the x-axis (C-scan). The second image (at γ2) was registered with the first image (at γ1) to calculate both M(θ, ϕ, β) and M′(θ′, ϕ′, β′) in the laboratory coordinates (Eqs. (5) and (6)). For this reason, the second image volume including the M′ were rotated back by γ1˗ γ2 around the x-axis. Due to the optical refraction effect at the tissue surface, the incident light deviated from the initial incident direction (the z-axis) which was defined by the A-scan direction of the imaging system. Such an effect changed the orientation of the “evaluation” plane and distorted the original 3D rectangular image grid (i, j, k) and was corrected before image registration [17]. The sample surface boundary was first identified using an intensity-based thresholding method as in our previous studies [13,14]. A 5 × 5 pixel median filter was applied to remove noisy pixels on the resulting 2D surface. The surface normal vector was then calculated using the “surfnorm” function in Matlab. The light direction (θ, ϕ) inside the sample after refraction was calculated using Eq. (3). The actual pixel coordinates (i′, j′, k′) for the original pixel (i, j, k) were calculated as:
A new rectangular image matrix was then constructed using 3D cubic spline interpolation based on the actual pixel coordinates (i′, j′, k′) obtained in Eq. (8). And the M vector was measured for each pixel using Eq. (5). This refraction correction procedure was applied to both image volumes acquired at two different angles in the dual-angle imaging procedure. Then the corrected image volume obtained at the second imaging angle was rotated back to register with that obtained at the first imaging angle. The 3D fiber orientation was then calculated for each registered pixel using Eq. (6). All image processing was implemented in Matlab. The 3DSlicer software (www.slicer.org) was used for visualizing the 3D OPT results.3. Results
The dual-angle OPT method was first verified by imaging the 3D fiber orientation in the mouse extensor digitorum longus (EDL) muscle that was placed at various known orientations. Once the 3D fiber direction was measured at one sample position, it can be used to predict fiber direction at other sample positions. Therefore, the results obtained in dual-angle OPT can be validated by comparing the expected values with the measured values at various sample positions. The accuracy and capability of dual-angle OPT was further investigated and analyzed by imaging the 3D fiber structure in the mouse tibialis anterior (TA) muscle and the bovine articular cartilage.
3.1 Extensor digitorum longus (EDL) muscle
Figure 2(a) illustrates the geometry of the EDL muscle in the dual-angle OPT test. The muscle sample was excised from a C57BL/6 mouse and fixed in 4% paraformaldehyde. The muscle sample was mounted on a two-axis rotation stage carefully aligned with the scanning directions to provide accurate rotations around the x-axis (C-scan) and y-axis (B-scan). To implement the dual-angle scanning, the EDL muscle was first scanned at γ1 = −20° (rotation around the x-axis); then imaged again at γ2 =+ 20°after rotating the sample by 40° around the x-axis. Figure 2(b) shows the conventional single-scan OPT of the EDL muscle at two different angles (γ1 and γ2) which were later used for dual-angle reconstruction. The fiber orientation measured in the single-scan OPT was the 2D projected angle of the 3D fiber in the evaluation plane (Fig. 1). As described in previous studies [14, 15], such projected tractography was constructed within the equi-transmural layer from the tissue surface.
Fig. 2 Single-scan OPT and registration of two image volumes acquired in dual-angle OPT. (a) A schematic illustration of the geometry of the dual-angle OPT imaging. The EDL sample was rotated in the laboratory coordinates by γ1 and γ2 around the x-axis to complete the dual angle imaging. (b) The single-scan OPT of the EDL muscle obtained at γ1 = −20° and γ2 =+ 20°. (c) The surface profile of the EDL muscle imaged at γ1 = −20°. (c) The surface profiles of the same EDL muscle imaged at γ2 =+ 20° and after rotating back by −40° around the x-axis to be registered with the results in (c). The unit in the contour plots (c) and (d) represents the height of the sample surface (in “pixels” along the z-axis).
Before dual-angle image construction, the registration of the two image volumes measured at γ1 and γ2 was examined by comparing their surface profiles. As expected, the “contour” profiles of the sample surface obtained at γ1 = −20° (Fig. 2(c)) and γ2 =+ 20° (Fig. 2(d)) were noticeably different due to their different image positions. To register the sample surface measured at γ2 =+ 20° with that obtained at γ1 = −20°, the image data acquired at γ2 = 20° was rotated back by 40° by multiplying Rx(−40°) to the 3D coordinates of each image pixel, where Rx is the standard 3D rotational matrix around the x-axis. The resulting surface (Fig. 2(d), right) appeared very similar to the one shown in Fig. 2(c). Due to the alignment errors of the rotational stages, the optimal rotation angle was slightly adjusted within a small range (<2°) to minimize the pixel-wise error between the two surface profiles. Quantitatively, 95% of the sample surface had < 5-pixel difference after registration.
During the validation process, the above dual-angle OPT procedure was performed multiple times while placing the EDL muscle at six different positions from α = −30° to + 20° in a step of 10° around a “validation” axis V in the y-z plane (Fig. 3(a)). The rotational axis V had a polar angle θv = 90°˗ γ1 with the z-axis. Two example 3D image volumes obtained at α = −30° and + 20° were illustrated in Fig. 3(a) using different colors. These six positions induced considerable changes of fiber orientation in the three orthogonal projection planes (x-z, z-y, and x-y for AC, AB, and BC, respectively). The coordinates and axis V shown in Fig. 3(a) were actually located inside the tissue during the experiment (see Fig. 2(a)).
Fig. 3 3D dual-angle OPT and registration of multiple image volumes obtained at different 3D positions. (a) A schematic illustration of rotating the EDL muscle from α = −30° to α =+ 20° around the validation axis (V). (b) The 3D tractography of the EDL muscle placed at α = 0°, α =+ 20°, and α = −30° obtained using the dual-angle OPT at each position. (c) The surface profiles of the sample located at α = 0°. The images in (d) and (e) show the corresponding surface profile measured at α =+ 20° and α = −30°, and after rotating the image volumes back by −20° and + 30° around the validation axis to be registered with (c). The unit in the contour plots (c), (d), and (e) represents the height of the sample surface (in “pixels” along the z-axis).
At each sample position α, 3D fiber orientation was constructed using dual-angle OPT obtained at γ1 = −20° and γ2 =+ 20° following the procedure described in the methods section. Figure 3(b) shows example 3D tractography when the EDL muscle was placed at α = 0°, + 20°, and −30°. Different from the 2D projected angle shown in Fig. 2(b), the tractography in Fig. 3(b) represented the true 3D fiber orientation in space. These images constructed at different positions (α) were also registered to each other using the same procedure applied in registering the dual-angle images.
Figures 3(c) to 3(e) show the tissue surface imaged at α = 0°, + 20°, and −30° respectively. All these surface contours were from the first position of the dual-angle imaging (γ1 = −20°). These surface contour patterns showed a significant difference because the EDL muscle was imaged from different angles. However, after rotating the 3D volume back around the axis V (Fig. 3(a)), the resulting surface contours in Figs. 3(d) and 3(e) appeared quite similar to the sample surface imaged at α = 0° (Fig. 3(c)). In all cases, 92% of the surface area had < 5-pixel difference after rotation, which confirmed the goodness of the image registration.
For the purpose of validation, the directly measured 3D fiber angle at a sample position α was compared with the “expected” angles calculated using 3D fiber angle measured at α = 0°. The “expected” 3D fiber direction Fe(α) = [fx(α), fy(α), fz(α)]T at sample position α can be calculated from the 3D fiber direction F(0) measured at sample position α = 0:
where Rx and Ry are standard rotational matrices around the x- and y-axis, respectively. To facilitate the comparison among 3D angles, we compared the 2D projected angles of the 3D direction in three orthogonal projection planes calculated in Eq. (7).In practice, the 3D angle measured at α = 0° had experimental errors, which may lead to inaccurate angles calculated at other sample positions. Therefore a “genetic algorithm” based optimization procedure was used to find the best 3D fiber direction at α = 0°. The original fiber orientation obtained in dual-angle OPT was used as the initial value. True fiber orientation was obtained in the optimization to minimize an error function defined as:
where stands for the expected fiber angle in the projection plane p calculated by Eq. (9) and denotes the measured projection angle from dual-angel OPT at the same projection plane. The optimized 3D fiber angle at sample position of α = 0° was then used to calculate the expected projection angels at other sample positions.Figures 4(a)-4(c) show an example of optimization of the fiber orientation measured in a 48 × 48 × 48 µm3 region-of-interest (ROI). For this ROI, the polar and azimuthal angles measured from dual-angle OPT at α = 0° were 73.2° and 28.7°, respectively; whereas the corresponding optimized angles were 70.8° and 28.8°. Since the difference in the azimuthal angle was negligible, the overall error (Eq. (10)) was slightly greater in the AC-plane (0.73° ± 2.13°) and AB-plane (−0.25° ± 2.24°) than the BC-plane (0.22° ± 0.98°). Rotating the EDL muscle from α = −30° to α = 20° changed the fiber orientation from 6° to 35° in the AC-plane, −6° to −62° in the AB plane, and 43° to 20° in the BC-plane.
Fig. 4 Example optimization results. The symbols represented the measured projection angles for a region-of-interest (ROI) in the AC (a), AB (b), and BC (c) planes at different sample positions. The solid line indicated the optimization results (best-fitted) based on Eq. (10).
Figure 5 shows the validation results obtained from 50 randomly selected ROIs (size 48 × 48 × 48 µm3) in a 0.3 × 0.3 × 1.0 mm3 (A × B × C) volume of the EDL muscle. Since each ROI was measured at six different sample positions from α = −30° to + 20°, there were 300 data points in the correlation analysis between the measured and “expected” projection angles in the AB-, AC- and BC-plane. In all three projection planes, the “expected” angles were highly correlated with the actual measured angles (Fig. 5). All three linear regression results showed a close to 1.0 slope and ≤0.51° y-intercept. The coefficients of determination (R2) were 0.89, 0.95, and 0.94 in the AC-, AB-, and BC-plane, respectively. These results suggested that the experimentally measured 3D fiber orientation in dual-angle OPT closely matched with the “expected” results. Among all data points, the mean and standard deviation of the difference between the measured and expected fiber angle were 0.00° ± 3.30°, 0.00° ± 4.02°, and 0.14° ± 1.96° in the AC, AB, and BC planes, respectively.
Fig. 5 The correlation between measured and expected projection angles in fifty randomly selected ROIs (48 × 48 × 48 µm3) from the EDL muscle. The projection angles were calculated in (a) AC-plane, (b) AB-plane, and (c) BC-plane. The regression equation and the coefficient of determination (R2) were also shown. The distribution of (d) coefficient of determination (R2), (e) slope, and (f) y-intercept of the linear regression analyses between measured and expected fiber angles obtained in ten tests of using 50 random ROIs in the EDL muscle.
Similar results were obtained when different random sets of ROIs were used in the validation. We repeated ten times the aforementioned procedure of selecting 50 ROIs randomly in the EDL muscle. The mean and standard deviation of the coefficients of determination (R2) were 0.88 ± 0.01, 0.95 ± 0.01, and 0.94 ± 0.01 in the AC-, AB-, and BC-plane, respectively (Fig. 5(d)). The corresponding slopes of the linear regression were 0.99 ± 0.02, 0.98 ± 0.01, and 0.98 ± 0.02 (Fig. 5(e)); and the y-intercept values were 0.10 ± 0.46, −0.72 ± 0.44, and 0.19 ± 0.77 (Fig. 5(f)) in the AC-, AB-, and BC-plane respectively. These results indicated a high consistence between the results of dual-angle OPT measurement and the expected values and thus validated dual-angle OPT for accurately imaging 3D fiber orientation.
3.2 Tibialis anterior (TA) muscle
In this section, we applied dual-angle OPT to image the 3D fiber organization in the mouse TA muscle. The muscle was excised from a C57BL/6 mouse and fixed in 4% paraformaldehyde. The TA sample was imaged at γ1 = −25° and γ2 =+ 25° in dual-angle OPT.
Figure 6(a) shows the 3D intensity images of the TA muscle imaged at γ2 =+ 25°. Figure 6(b) shows the corresponding 3D tractography constructed in dual-angle OPT. The 3D tractography indicated that myofibers were mostly aligned along the long axis of the TA muscle. Figure 6(c) demonstrates the capability of 3D OPT in revealing the detailed fiber architecture by selectively revealing two fiber bundles at two separate locations inside the muscle. As expected, these two fiber bundles were clearly organized along the long axis of the muscle toward the apex. Figure 6(d) further demonstrates that the details of the fiber bundles can be studied by projecting the 3D fiber bundles to 2D projection planes (the AB-, AC-, and BC-plane).
Fig. 6 Dual-angle OPT of a piece of the mouse tibialis anterior (TA) muscle. (a) The 3D intensity image of the muscle. (b) The 3D tractography obtained in dual-angle OPT. (c) Detailed fiber organization of two separate fiber bundles inside the muscle. (d) The 2D tractography by projecting the 3D fibers shown in (c) to the AB, AC, and BC planes.
Because of its large size, the TA muscle fibers were still discernable in the intensity images. Therefore the intensity images in this case can serve as another validation standard to evaluate the dual-angle OPT. As shown in Fig. 6, the 3D fiber orientation obtained in dual-angle OPT can be projected to any 2D planes in space. For the purpose of this validation, we chose to compare the fiber orientations in the en face (BC) plane. The 3D fiber vectors were projected onto the en face image plane of the sample placed at the second dual-angle position (γ2 =+ 25°). The fiber orientations at the same en face plane were also calculated from the corresponding intensity images as described previously [17]. Briefly, the pixel-wise 2D gradient was calculated using a 3 × 3 Sobel window. The distribution of the gradient direction was fitted using the Von Mises distribution within a 200 × 200 µm2 evaluation window. The orientation corresponding to the peak of the fitter distribution was assigned as the primary fiber orientation at the center of the evaluation window.
Figure 7 shows the comparison of the projection of 3D OPT results with the intensity tractography in the en face plane of the sample imaged at γ2 at different depths up to 800 μm. Overall, the dual-angle OPT projection was consistent with intensity based tractography. Myofibers converged along the long muscle axis toward the apex at all depths from the intensity images, which was also evident in the dual-angle OPT results. However, there were also some major discrepancies between these two measurements. The intensity-based tractography appeared to be much noisier, especially at the upper-right part, than the OPT results. The edge-detection based intensity tractography is sensitive to any intensity variations. Such variations may not be related to muscle fibers, but can induce artificial fiber orientation changes in the tractography. At depths greater than ~600 µm, the striated fiber patterns can no longer be clearly resolved in the intensity image. Therefore, the intensity tractography showed larger variations; whereas dual-angle OPT still showed relatively a smooth trend of the fibers toward the apex. When overlaying the OPT results on top of the intensity images, we verified that dual-angle OPT provided a better match with the fiber trajectories than the intensity based method.
Fig. 7 A comparison of the fiber organization obtained in dual-angle OPT with that obtained by analyzing OCT intensity profiles in the en face plane at depths from 200 µm to 800 µm. The first row shows the en face intensity obtained in the muscle sample place at γ2 =+ 25° in dual-angle OPT. The second row shows the tractography calculated from the intensity images. The third row shows the 2D tractography calculated by projecting the 3D fiber orientation obtained in dual-angle OPT to the corresponding en face plane (Eq. (7)).
3.3 Articular cartilage
In this section, dual-angle OPT was applied to image the collagen fiber organization in the bovine articular cartilage. The cartilage sample was excised from the middle phalanx and fixed in 4% paraformaldehyde. This sample was imaged from the cartilage surface (referred to as the “top-scan”) using the dual-angle procedures at γ1 = −25° first and then at γ2 =+ 25°.
No fibers can be observed from the 3D intensity image of the cartilage (Fig. 8(a)). However, the constructed 3D fiber tractography from dual-angle OPT (Fig. 8(b)) clearly revealed the unique arcade-like fiber organization in the cartilage [26]. In the “superficial” zone, fibers were oriented relatively parallel to the sample surface. Then the fibers arched in the “transitional” zone to be eventually perpendicular toward the boundary between non-calcified and calcified cartilage in the deeper “radial” zones. Moreover, the 3D tractography correctly revealed the “brushing” direction, i.e. all fibers were bended toward the same direction [22].
Fig. 8 3D visualization of a piece of the bovine cartilage in (a) intensity, (b) 3D fiber tractography from dual-angle measurements, and (c) the 2D tractography by projecting the 3D fibers shown in (b) to the AB-, AC-, and BC-plane as labeled in (a).
Similar to the example of the TA muscle (Fig. 7), the 3D OPT data set can be explored programmingly to provide a detailed view of the fiber architecture. Due to the unique “arcade” 3D fiber structure in the cartilage, the fiber orientation would appear differently when viewed from different angles or in cutting planes made at different orientations. Figure 8(c) demonstrated such an effect by showing the fiber organization when projecting the 3D orientation to three orthogonal 2D cutting planes (Eq. (7). This result underscored the importance of considering the cutting directions when evaluating the histology and scanning electron microscopy (SEM) images of cartilage samples [26].
To verify the obtained 3D fiber architecture, one side of the cartilage (indicated in Fig. 8(a)) was directly imaged using OPT. The intensity image (Fig. 9(a)) showed the boundary between the non-calcified and calcified cartilage, but no fiber structure can be observed. The resulted 2D tractography (Fig. 9(a)) clearly uncovered the arcade-like fiber structure that was conventionally used to identify three different zones in the cartilage [26–28]. The fiber appeared to be randomized in the calcified region. As a comparison, the 3D fiber orientation data measured at the same location using dual-angle OPT were projected on the same side-scan plane. Although this plane was parallel to the incident light direction or perpendicular to the OPT evaluation plane (Fig. 1), the projected 2D tractography (Fig. 9(b)) appeared qualitatively in good agreement with the directly measured “side-scan” result (Fig. 9(a)).
Fig. 9 (a) The en face intensity and tractography images obtained from direct OPT imaging of the side of the cartilage sample (the “side-scan” as labeled in Fig. 8). (b) The tractography by projecting the 3D fiber orientation obtained in dual-angle OPT to the same tissue side as imaged in (a). NC: non-calcified cartilage; CB: calcified cartilage/bone.
We further quantitatively compared the “side-scan” results with the 3D projection. Figure 10(a) shows the depth-dependent fiber orientation obtained by averaging all line-scans within a ROI (marked with a red box in Fig. 9(a)). The “transitional” segment where the fiber orientation changed from parallel to perpendicular can be well fitted (adjusted R2 = 0.999) using a hyperbolic tangent function as reported in a previous study using polarization light microscopy (PLM) [27]. The fitted equation, y = 42.06 × tanh[(x-121.8)/54.47] + 46.13, was plotted as a blue line in Fig. 10(a). The first derivative of the fitting equation was shown as the dark-green curve. The boundaries of the transitional zone were determined as the boundary positions of the full-width-half-maximum (FWHM) of the first derivative curve [28]. This procedure resulted in ~74 µm as the boundary between the superficial and transitional zone, and ~170 µm as the boundary between the transitional and the radial zone.
Fig. 10 (a) The fiber orientation over the depth obtained from the side-scan image by averaging all B-scan inside a 0.5 mm wide ROI (marked in red box in Fig. 9(a)). (b) The orientation change over depth obtained from the projected 2D fiber orientation (Fig. 9(b)). The blue lines in (a) and (b) are fitting results using a hyperbolic tangent function. The green curves are the first order derivative of the curve fitting. The labels “ST”, “TR”, and “CB” indicate the boundaries between the superficial and transitional zone, the transitional and radial zone, and cartilage and bone, respectively. Error bars indicate standard deviation.
The above procedure was repeated using the projected angle θp calculated from the 3D fiber orientation using dual-angle OPT (Fig. 9(b)). We used n = 1.5 as the optical refractive index of the cartilage as reported in a previous study [29]. The orientation profile can also be well fitted (adjusted R2 = 0.989) using a hyperbolic tangent function: y = 34.69 × tanh[(x-113.7)/40.14] + 44.7. Using the same criteria described in the previous paragraph, the boundary between the superficial zone and the transitional zone was located at ~79 µm which was consistent with that determined from direct side-scan. The boundary between the transitional zone and the radial zone was at ~150 µm which was slightly smaller than that determined from the side-scan. It was noticed that the fiber orientation had larger fluctuation at the junction between the transitional zone and the radial zone, which may have caused larger measurement errors in both cases. Similar observation has been reported previously using PLM [27], which was attributed to fiber variation in the transitional zone.
The directly measured side-scan results (Fig. 10(a)) indicated that the orientation started to decrease significantly at ~660 µm. Similarly, the projected angle from 3D measurements also showed a significant decrease in fiber orientation at ~656 µm. This location was likely the boundary between the non-calcified cartilage and calcified cartilage or bone. The fiber orientation also had large variations within this region, which may explain the difference between the two measurements within this region. At locations inside the bone region, the side-scan results showed better imaging depth than the 3D projection results.
4. Discussions
Imaging 3D fiber organization is important to understand mechanical properties and normal physiological functions in fibrous tissues. In this study, 3D fiber orientation was successfully imaged by using a dual-angle OPT procedure. The proposed method utilizes tissue optical polarization properties to calculate actual 3D fiber orientations. The accuracy and performance of dual-angle OPT were investigated by testing three different biological tissues. These studies indicated that dual-angle OPT can reliably acquire the full 3D images of microscopic fiber organization in complicated biological tissues.
The EDL muscle was used in the validation test due to its uniform and well-organized fiber structure. The experiment was carefully designed to ensure the rotation axis and the amount of rotation can induce sufficient changes in fiber orientation. The quantitative comparison between the expected and the experimental values showed consistent accuracy of dual-angle OPT. The two images acquired in the dual-angle test were accurately registered at each sample position. In addition, the sample images acquired at multiple different rotational positions were also registered to the reference location at α = 0°. However, it is challenging to achieve a perfect registration in all 12 images used in the validation. We noticed that a couple of distant data points had relatively large errors in the regression results (Fig. 5). A close examination revealed that these points were part of the dual-angle imaging at the two largest rotational positions at α =+ 20° and α = −30°. This error was likely caused by the large surface inclination at these sample positions. A larger light incident angle induced bigger beam deviation, which was more sensitive to slight errors in image registration.
Although the measured and expected fiber orientations were highly correlated in all three orthogonal projection planes, the correlation in the AC projection plane was marginally poorer than in the other two planes (Fig. 5). We believe this can be attributed to two factors after analyzing the optimization results obtained in all ROIs, First, we found that the measured polar angle of the 3D fibers had a very small (3.15° ± 1.75°) but consistent difference from the optimized results; whereas the corresponding azimuthal angle had better agreement with the optimized values (1.60° ± 1.57°). As shown in Eq. (7), a larger variation in the polar angle θf affected the z-axis related projection planes (AB- and AC-plane) more than the BC-plane. Such a slightly bigger variation in the polar angle was likely due to an imperfect alignment of the rotational axis. Second, it is known that a similar amount of error in the regression model may lead to a smaller R2 if the measurements had a smaller range. Fiber orientations changed much greater in the AB-plane (~80°) than in the projection planes associated with the x-axis (~40° changes in the AC- and BC-plane). Therefore the measurements in the AC-plane were affected by both factors as discussed above, which resulted in a slightly worse correlation.
The TA muscle was used as another validation because the muscle fibers were discernable in the intensity images as well. This is because the muscle fibers were bigger in the TA muscle than in the EDL muscle as revealed in a histology study [30]. In general, the fiber orientation obtained in dual-angle OPT was consistent with that obtained using intensity based image processing. Such a direct comparison also revealed some important advantages in the OPT method. In order to resolve fibers using intensity based contrast, the image resolution should be sufficiently high to detect intensity changes around the fibers or fiber bundles. Then image processing methods were commonly applied within an evaluation window to determine the primary orientation inside the window. The size of the evaluation window sometimes reached a few hundred-micrometers [8]. Therefore, the actually spatial resolution was compromised in computing the fiber orientation. In addition, the intensity based image processing may be affected by any intensity changes that are not related to fibers. For example, OCT intensity images may contain the “banding” artifacts in a birefringent tissue caused by the residual polarization effects in the imaging system [31]. Such non-fiber related structural changes may cause incorrect fiber calculations. This issue became more severe at large depths when intensity contrast deteriorated as shown in Fig. 7.
Intensity-based tractography was commonly applied to the en face images [6,8], and thus it can obtain the projected fiber orientation in the en face plane. A computational approach was previously explored to calculate the 3D fiber orientation from the measured 2D orientation in the en face plane [9]. It relied on the assumption that the fibers were parallel to the sample surface, which is not generally applicable for most fibrous tissues including heart, cartilage, and neural tissues. A structure tensor based algorithm was recently reported to measure 3D fiber orientation by analyzing the 3D intensity variation [10]. However, this approach has the same drawback as in any intensity based methods as described above.
As a comparison, OPT can maintain the native spatial resolution of the imaging system. It calculates fiber orientation using the local optic axis on a single-pixel basis. The optic axis calculation in OPT is not affected by structural variations in the intensity image. OPT does not require a spatial resolution that is high enough to resolve individual fibers or fiber bundles. In fact, the intensity images failed to reveal any fibers in both the EDL muscle and bovine cartilage. However, OPT still successfully imaged the local fiber orientation in these samples, which represented the averaged fiber orientation at each image pixel. The dual-angle OPT procedure described in this work relied upon basic geometrical transformations and did not need any presumptuous information regarding the fiber structure. The results shown in Fig. 7 supported that OPT worked better than intensity-based method especially at larger imaging depths.
The articular cartilage results best demonstrated the advantage of the proposed dual-angle OPT. Articular cartilage is structurally reinforced by a network of collagen fibers which are organized in a special “arcade” architecture [22, 26]. A single-scan OPT can still map the fiber orientation in the cartilage similar to imaging applications in other tissue samples [14–16]. However, single-scan OPT measures the projected fiber angle in the 2D evaluation plane as illustrated in Fig. 1. These angles measured in single-scan OPT may correctly represent the fiber orientation in the superficial cartilage where collagen fibers are tangential to the sample surface. The fibers are almost perpendicular to the surface or aligned with the incident light in the deep radial zone of the cartilage. Therefore the measured 2D projection angles in the radial zone do not represent the complete picture of the true 3D fiber orientation there.
Dual-angle OPT clearly revealed the well-known arcade-like fiber organization in the cartilage. The obtained 3D fiber orientation can be projected to the 2D side plane of the cartilage and compared with the directly measured fiber orientation. Despite the challenge of achieving a perfect image registration, the projected fiber orientation is consistent with the directly measured fiber orientation. The unique fiber organization in the cartilage plays a critical role in maintaining its proper mechanical properties. Changes in collagen fiber structure including the orientation are signs of early cartilage degeneration [32]. Therefore imaging fiber organization in cartilage may be useful for detecting early cartilage diseases. Polarization light microscopy (PLM) [27, 28] is currently the standard way for imaging fiber orientation in cartilage. Unfortunately, PLM can only image thin cartilage slices due to its limited imaging depth. Being able to map the true 3D fiber structure in cartilage, dual-angle OPT provides a better nondestructive alternative to PLM and may find many valuable applications in cartilage studies.
The geometric transformations used in dual-angle OPT rely on the local fiber orientation measured in each individual OPT as well as the orientation and location of the corresponding “evaluation plane” (Fig. 1). It is important to note that the “cumulative” optic axis measured in conventional PSOCT studies does not represent the correct fiber orientation in tissue [11] and cannot be used for constructing the 3D orientation. The surface refraction changes the orientation and location of the evaluation plane (Fig. 1). Although its effect may not be significant at small incident angles or in relatively homogeneous samples, correcting its effect can ensure a proper image registration for dual-angle calculation in general situation. The impact of refraction correction was observed in the validation results using EDL muscle due to the large fiber orientation changes induced by sample rotation. Following the exact same procedures, but using a refractive index of 1.0, we found that the correlation results shown in Fig. 5 were negatively affected especially in the BC projection where the correlation R2 decreased to 0.79, the fitting slope decreased to 0.74, and the intercept increased to 8.99°.
The surface refraction may also introduce a polarization change because of the different transmission in the s- and p-polarized components. As an example, such a transmission difference is ~3.3% at a 30° incident angle and 1.5 sample refractive index based on Fresnel’s equation. Because this small difference introduced a <1° change in polarization orientation, it was not considered in the current study. The current refraction correction method may require a relatively smooth sample surface. For tissues with rough surface, refractive index matching solution might be used to circumvent the need of refraction correction.
It may be beneficial to use a small incident angle as a large incident angle may also affect the effective imaging depth. For example, a 1-mm imaging depth may be reduced to 866 µm with a 30° refractive angle due to oblique optical path. Dual-angle OPT only calculates the fiber orientation in the common image volumes acquired at two different sample positions. Therefore its effective image depth is limited by the bigger refractive angle between the two volumes. In this study, the imaging depth was greater than 800 µm in the TA muscle and ~700 µm in the cartilage sample when ± 25° incident angles were used in the dual-angle measurements. Limiting the incident angle can also reduce the aforementioned refraction induced polarization changes. Nevertheless, using a light source at a longer wavelength (e.g. 1.3 μm) may achieve greater imaging depths in each single-scan OPT and thus can extend the overall effective imaging depth in dual-angle OPT.
Finally, the need of two images from the same sample acquired at different angles might be challenging for in vivo applications because sample movement may reduce the size of the common volume used in dual-angle reconstruction. Therefore, innovative system designs that can achieve simultaneous dual-angle scanning might be desirable for in vivo imaging.
5. Conclusion
We proposed a method to nondestructively image 3D fiber orientation in a tissue. This method utilized the recently developed OPT method at two different imaging angles in relation to the sample. Each OPT measured the projected 2D orientation within a different evaluation plane. The 3D fiber orientation was then constructed using geometrical transformation based on the two measurements. We described an example implementation of this method and conducted several tests to verify the obtained 3D angles in several different biological tissues. In particular, dual-angle OPT successfully imaged the unique depth-varying fiber orientation in cartilage. As an enabling tool, dual-angle OPT provides a promising platform for studying the important 3D fiber architecture in many fibrous tissues.
References and links
1. B. Taccardi, E. Macchi, R. L. Lux, P. R. Ershler, S. Spaggiari, S. Baruffi, and Y. Vyhmeister, “Effect of myocardial fiber direction on epicardial potentials,” Circulation 90(6), 3076–3090 (1994). [CrossRef] [PubMed]
2. G. Buckberg, J. I. E. Hoffman, A. Mahajan, S. Saleh, and C. Coghlan, “Cardiac mechanics revisited: the relationship of cardiac architecture to ventricular function,” Circulation 118(24), 2571–2587 (2008). [CrossRef] [PubMed]
3. C. C. Van Donkelaar, L. J. G. Kretzers, P. H. M. Bovendeerd, L. M. A. Lataster, K. Nicolay, J. D. Janssen, and M. R. Drost, “Diffusion tensor imaging in biomechanical studies of skeletal muscle function,” J. Anat. 194(1), 79–88 (1999). [CrossRef] [PubMed]
4. D. D. Streeter Jr, H. M. Spotnitz, D. P. Patel, J. Ross Jr, and E. H. Sonnenblick, “Fiber Orientation in the Canine Left Ventricle During Diastole and Systole,” Circ. Res. 24(3), 339–347 (1969). [CrossRef] [PubMed]
5. P. J. Basser, J. Mattiello, and D. LeBihan, “MR diffusion tensor spectroscopy and imaging,” Biophys. J. 66(1), 259–267 (1994). [CrossRef] [PubMed]
6. C. P. Fleming, C. M. Ripplinger, B. Webb, I. R. Efimov, and A. M. Rollins, “Quantification of cardiac fiber orientation using optical coherence tomography,” J. Biomed. Opt. 13(3), 030505 (2008). [CrossRef] [PubMed]
7. C. M. Ambrosi, V. V. Fedorov, R. B. Schuessler, A. M. Rollins, and I. R. Efimov, “Quantification of fiber orientation in the canine atrial pacemaker complex using optical coherence tomography,” J. Biomed. Opt. 17(7), 071309 (2012). [CrossRef] [PubMed]
8. C. J. Goergen, H. Radhakrishnan, S. Sakadžić, E. T. Mandeville, E. H. Lo, D. E. Sosnovik, and V. J. Srinivasan, “Optical coherence tractography using intrinsic contrast,” Opt. Lett. 37(18), 3882–3884 (2012). [CrossRef] [PubMed]
9. Y. Gan and C. P. Fleming, “Extracting three-dimensional orientation and tractography of myofibers using optical coherence tomography,” Biomed. Opt. Express 4(10), 2150–2165 (2013). [CrossRef] [PubMed]
10. H. Wang, C. Lenglet, and T. Akkin, “Structure tensor analysis of serial optical coherence scanner images for mapping fiber orientations and tractography in the brain,” J. Biomed. Opt. 20(3), 036003 (2015). [CrossRef] [PubMed]
11. C. Fan and G. Yao, “Mapping local retardance in birefringent samples using polarization sensitive optical coherence tomography,” Opt. Lett. 37(9), 1415–1417 (2012). [CrossRef] [PubMed]
12. C. Fan and G. Yao, “Mapping local optical axis in birefringent samples using polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 17(11), 110501 (2012). [CrossRef] [PubMed]
13. C. Fan and G. Yao, “Imaging myocardial fiber orientation using polarization sensitive optical coherence tomography,” Biomed. Opt. Express 4(3), 460–465 (2013). [CrossRef] [PubMed]
14. Y. Wang and G. Yao, “Optical tractography of the mouse heart using polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 4(11), 2540–2545 (2013). [CrossRef] [PubMed]
15. Y. Wang, K. Zhang, N. B. Wasala, D. Duan, and G. Yao, “Optical polarization tractography revealed significant fiber disarray in skeletal muscles of a mouse model for Duchenne muscular dystrophy,” Biomed. Opt. Express 6(2), 347–352 (2015). [CrossRef] [PubMed]
16. L. Azinfar, M. Ravanfar, Y. Wang, K. Zhang, D. Duan, and G. Yao, “High resolution imaging of the fibrous microstructure in bovine common carotid artery using optical polarization tractography,” J Biophotonics (online 11 DEC 2015).
17. Y. Wang, K. Zhang, N. B. Wasala, X. Yao, D. Duan, and G. Yao, “Histology validation of mapping depth-resolved cardiac fiber orientation in fresh mouse heart using optical polarization tractography,” Biomed. Opt. Express 5(8), 2843–2855 (2014). [CrossRef] [PubMed]
18. M. A. Wallenburg, M. F. G. Wood, N. Ghosh, and I. A. Vitkin, “Polarimetry-based method to extract geometry-independent metrics of tissue anisotropy,” Opt. Lett. 35(15), 2570–2572 (2010). [CrossRef] [PubMed]
19. F. Fanjul-Vélez and J. L. Arce-Diego, “Polarimetry of birefringent biological tissues with arbitrary fibril orientation and variable incidence angle,” Opt. Lett. 35(8), 1163–1165 (2010). [CrossRef] [PubMed]
20. N. Ugryumova, S. V. Gangnus, and S. J. Matcher, “Three-dimensional optic axis determination using variable-incidence-angle polarization-optical coherence tomography,” Opt. Lett. 31(15), 2305–2307 (2006). [CrossRef] [PubMed]
21. D. K. Kasaragod, Z. Lu, J. Jacobs, and S. J. Matcher, “Experimental validation of an extended Jones matrix calculus model to study the 3D structural orientation of the collagen fibers in articular cartilage using polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 3(3), 378–387 (2012). [CrossRef] [PubMed]
22. Z. Lu, D. Kasaragod, and S. J. Matcher, “Conical scan polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 5(3), 752–762 (2014). [CrossRef] [PubMed]
23. C. J. Liu, A. J. Black, H. Wang, and T. Akkin, “Quantifying three-dimensional optic axis using polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 21(7), 070501 (2016). [CrossRef] [PubMed]
24. C. Fan and G. Yao, “Single camera spectral domain polarization-sensitive optical coherence tomography using offset B-scan modulation,” Opt. Express 18(7), 7281–7287 (2010). [CrossRef] [PubMed]
25. C. Fan and G. Yao, “Full-range spectral domain Jones matrix optical coherence tomography using a single spectral camera,” Opt. Express 20(20), 22360–22371 (2012). [CrossRef] [PubMed]
26. J. M. Clark, “The organisation of collagen fibrils in the superficial zones of articular cartilage,” J. Anat. 171, 117–130 (1990). [PubMed]
27. Y. Xia, J. B. Moody, N. Burton-Wurster, and G. Lust, “Quantitative in situ correlation between microscopic MRI and polarized light microscopy studies of articular cartilage,” Osteoarthritis Cartilage 9(5), 393–406 (2001). [CrossRef] [PubMed]
28. J. H. Lee and Y. Xia, “Quantitative zonal differentiation of articular cartilage by microscopic magnetic resonance imaging, polarized light microscopy, and Fourier-transform infrared imaging,” Microsc. Res. Tech. 76(6), 625–632 (2013). [CrossRef] [PubMed]
29. J. M. Herrmann, C. Pitris, B. E. Bouma, S. A. Boppart, C. A. Jesser, D. L. Stamper, J. G. Fujimoto, and M. E. Brezinski, “High resolution imaging of normal and osteoarthritic cartilage with optical coherence tomography,” J. Rheumatol. 26(3), 627–635 (1999). [PubMed]
30. V. Augusto, C. R. Padovani, and G. E. R. Campos, “Skeletal muscle fiber types in C57BL6J mice,” Braz. J. Morphol. Sci. 21(2), 89–94 (2004).
31. X. Yang, D. Lorenser, R. A. McLaughlin, R. W. Kirk, M. Edmond, M. C. Simpson, M. D. Grounds, and D. D. Sampson, “Imaging deep skeletal muscle structure using a high-sensitivity ultrathin side-viewing optical coherence tomography needle probe,” Biomed. Opt. Express 5(1), 136–148 (2014). [CrossRef] [PubMed]
32. V. B. Shim, T. F. Besier, D. G. Lloyd, K. Mithraratne, and J. F. Fernandez, “The influence and biomechanical role of cartilage split line pattern on tibiofemoral cartilage stress distribution during the stance phase of gait,” Biomech. Model. Mechanobiol. 15(1), 195–204 (2016). [CrossRef] [PubMed]
