Computational multi-directional optical coherence tomography for visualizing the microstructural directionality of the tissue

: We demonstrate computational multi-directional optical coherence tomography (OCT) to assess the directional property of tissue microstructure. This method is the combination of phase-sensitive volumetric OCT imaging and post-signal processing. The latter comprises of two steps. The first step is an intensity-directional analysis, which determines the dominant en face fiber orientations. The second step is the phase-directional imaging, which reveals the sub-resolution depth-orientation of the microstructure. The feasibility of the method was tested by assessing muscle and tendon samples. Stripe patterns with several sizes were visualized in the phase-directional images. In order to interpret these images, the muscle and tendon structures were numerically modeled, and the phase-directional images were generated from the numerical model. The similarity of the experimental and numerical results suggested that the stripe patterns correspond to the muscle fiber bundle and its crimping.


Introduction
Tissue microstructure and its directionality are highly related to tissue functions, and hence its evaluation is important. For example, the microscopic architecture of myocardial fiber is related to the beating function [1] and electrical conduction [2] of the heart. In the brain, the nerve fiber pathway is associated with brain function [3]. In the eye, the organization of Henle's fiber layer would give diagnostic information on retinal diseases [4]. Collagen fiber architecture is linked to the mechanical properties of connective tissues such as tendon [5] and ligament [6].
Several methods have been used to examine these microstructures. Histology is the gold standard for evaluating microstructures, but sample preparation is time-consuming and causes microstructural shrinkage and distortion. Diffusion-tensor magnetic resonance imaging can visualize three-dimensional fiber organization, and has been applied to the heart [7,8] and brain [9], but the resolution is limited to sub-millimeter scale. Second-harmonic generation microscopy can directly image collagen fiber and has been applied to the heart [10], tendon [11], and ligament [12]. However, the imaging depth is limited to a few hundred microns.
Optical coherence tomography (OCT) [13] is a widely used modality providing noninvasive, high-resolution, and high-penetration imaging. OCT has been applied to microstructural analysis by combining with extensive image processing and with hardware extensions. The examples of the former is OCT intensity-based tractography. It uses intensity-gradient analysis [14,15] and/or Radon transform [16] to analyze the dominant fiber orientation in en face OCT. Three-dimensional (3-D) fiber visualization has also been demonstrated. Gan et al. extracted 3-D fiber structure of the heart [17], where they assumed the inclined angle of en face fiber was parallel to the surface. McLean et al. demonstrated 3-D tractography of the human uterus by applying intensity-gradient method to the en face plane and Radon transform to the cross-sectional OCT images [18]. Since these methods use only the intensity of the OCT images, its performance is limited by the resolution of the OCT system.
The example of the latter, the hardware extension of OCT, is polarization-sensitive OCT (PS-OCT) based tractogaphy. It uses polarization property of the tissue to identify the fiber orientation, and it has been used to visualize myocardial fibers [19] and nerve fibers of the brain [20]. Although this modality was successful, it requires PS-OCT hardware and is more complicated than standard non-polarization sensitive OCT.
Another example of hardware extension is multi-directional OCT (MD-OCT) [21][22][23]. It measures a sample from multiple directions and visualizes the directionality of the tissue by showing the intensity differences among the multi-directional OCT signals. This method has successfully visualized the directionality of tissue such as Henle's fiber layer that cannot be clearly measured with conventional OCT. However, this method is only for crude observation.
In MD-OCT, the tissue microstructure can be assessed in more details by increasing the number of measurement directions, and there are two ways to achieve this. One is increasing the number of measurements, but it will significantly increase the measurement time. The other is using multiple probe beams with multiple incident directions. Wartak et al. demonstrated simultaneous three-beam imaging in which the beams were incident on the sample from different directions [24]. Although it is not to assess the tissue microstructure, Blatter et al. and Klein et al. demonstrated two-beam [25] and four-beam imaging [26], respectively. However, these approaches requires complicated hardware, and also the power of each probe beam, and hence the sensitivity of each imaging channel, decreases as the number of probe beams increases. Hence, it is not realistic to further increase the number of probe beams, or similarly the number of measurement directions.
We previously demonstrated the method which realizes multi-directional imaging through a signal processing of complex OCT [27]. In this paper, we further enhanced this method, so denoted as computational multi-directional OCT. This method applies numerical sub-aperture masks to the spatial frequency spectrum of an en face OCT image. It realizes multi-directional imaging purely through post-measurement signal processing and provides microstructural information on the sample. Since the multi-directional images are generated after image acquisition, we can generate arbitrary numbers of directional OCT images from arbitrary directions without performing more measurements. In addition, this sub-aperture processing can be applied independently to the intensity and the phase of the en face OCT. The numerical sub-aperture masking of the OCT intensity gives the dominant en face fiber orientations, while that of the phase reveals the sub-resolution depth-orientation of the microstructure that is not visible in the OCT intensity image. In addition, since this method is a purely numerical method, it requires only a standard, i.e., non-polarization-sensitive OCT hardware. The feasibility of this method is demonstrated by assessing bovine Achilles tendon and chicken breast muscle samples.

Principles of computational multi-directional imaging
We propose a computational scheme of multi-directional imaging for analyzing tissue microstructures, which is a sequence of two types of numerical image processing. The first is "intensity-directional analysis," which processes the intensity of an en face OCT image and determines the dominant en face fiber orientation. The second is "phase-directional imaging," which processes the phase of the en face OCT image and reveals the depth-oriented property of the microstructure. The details of these methods are as follows.

Intensity-directional analysis
The intensity-directional analysis is for determining the dominant orientation of an en face OCT and is achieved by spatial frequency masking of the en face OCT as shown in Fig. 1. The intensity of the en face OCT [ Fig. 1(a)] is numerically two-dimensionally (2-D) fast Fourier transformed (FFTed) to its spatial frequency spectrum [ Fig. 1(b)]. In order to extract the spectral components of a specific microstructural orientation, we apply a numerical mask with a paired fan-shaped opening [ Fig. 1(c)], similar to Ref. [16], to the spectrum. The total squared intensity of the masked spectrum is computed to represent the signal energy of the specific orientation structure.
Each fan-shaped opening of the mask has a 20 • opening angle, and the two openings are symmetrically centered at the zero-frequency. The spatial frequency range of the mask is set from 5 to 50 mm −1 . The minimum frequency is chosen to meet the minimum spatial frequency of fibers in the OCT image. The maximum frequency is ten times the minimum frequency to have some extent of spatial frequency around the fibers. Example of an en face OCT image (a) and its spatial frequency spectrum (b), which is generated by numerical 2-D Fourier transform. A numerical mask (c) is applied to the spatial frequency spectrum, and the total squared intensity of the masked spectrum is computed. The total squared intensity is plotted as a function of the rotation angle of the mask (d). The dotted red line indicates the centroid of the peak.
To analyze the distribution of the fiber orientation, we compute the total squared intensity of the masked spectrum for several mask orientations. In the current specific implementation, the mask is rotated from 0 • to 179 • with a 1 • interval, so that 180 total squared intensities are computed. By considering the point symmetry of the mask, these 180 total squared intensities cover whole 360 • rotation. Here the 0 • orientation is defined as shown in Fig. 1(c). The total squared intensities of each orientation are plotted as a function of the rotation angle as shown in Fig. 1(d). In this paper, this plot is denoted as the "angle distribution of structural orientation." The dominant fiber orientation is finally defined as the orientation orthogonal to the peak angle. Note that the peak indicates the dominant orientation of the structural frequency, from which the dominant fiber orientation is determined to be 90 • from the peak.
In our implementation, the peak angle is defined as the centroid of the peak, which is computed as follows. In order to avoid angle wrapping, the maximum rotation angle is set to be the center of the angle range. The estimate of the centroid is computed within a ± 15 • range around the peak. It is empirically known that this ± 15 • range covers the spread of the peak when the opening angle of the mask is 20 • .

Phase-directional image processing
After applying the intensity-directional analysis, the OCT data is processed with a "phasedirectional imaging" algorithm to visualize the depth-oriented directionality of microstructures. Here a complex en face slice is extracted from the complex OCT volume. To process only the phase of the slice, its amplitude is set to unity. This en face phase slice is 2-D Fourier transformed to the spatial frequency domain. And a binary mask with a single fan-shaped opening is applied, exemplified in the insets of Figs. 2(b)-(e). The masked spectrum is inversely Fourier transformed, and its intensity is defined as an en face"raw phase-directional image." In our particular implementation, the center of a fan-shaped mask's arc is at the zero frequency, and the opening angle is 180 • . The opening covers the frequency range of 10 to 210 mm −1 , which occupies almost all of the spatial frequency range as the Nyquist frequency is 256 mm −1 . The minimum cut-off frequency of 10 mm −1 was chosen to exclude the non-structured component, i.e., the DC component, and hence to highlight the directionality of microstructures. The 180 • opening angle was chosen to maximize the opening area as keeping its directionality, and hence to best preserve the lateral resolution of the raw phase-directional image. This process is repeated for several mask orientations including parallel and perpendicular to the dominant fiber orientation identified in the intensity-directional analysis described in Section 2.1. Note that, for each parallel and perpendicular orientations, two fan-shaped masks in 180 • relation are used. Hence, four binary masks are used in total.
Examples of this process is shown in Fig. 2. Figure 2(a) is the intensity en face slice and whose phase is processed. Figures 2(b)-(e) show four raw phase-directional images created using four binary masks (insets) with four orientations.

Pseudo-color visualization
To visualize the difference between raw phase-directional images of two diagonal orientations, a pseudo-color phase-directional image is generated as schematically depicted in Fig. 3. This pseudo-color image is made from two diagonally oriented raw phase-directional images. Here one image is used for the red channel and the other is used for the green channel. The blue channel is the average of these two. Prior to this color composition, each raw phase-directional image is Gaussian blurred with a standard deviation of 2 pixels (3.9 µm) to reduce the speckle. And in order to further emphasize the directional difference, the color saturation of the pseudo-color phase-directional images is doubled. Finally, in order to visually correlate the phase-directional image and the OCT intensity image, the brightness (visibility) of the image is reset to be the en face OCT intensity.
As shown in the color bar on the right of Fig. 3, the pseudo-color image becomes purple or magenta if the red-color orientation is dominant, while the image becomes light green if the green-color orientation is dominant. The image becomes gray if there is no difference between two orientations. Fig. 3. Diagram of pseudo-color phase-directional image formation. Two raw phasedirectional images (a) and (b) are assigned to red (c) and green (d) channels respectively. The blue channel (e) is the average of these two raw phase-directional images. Finally, a pseudo-color phase-directional image is generated by merging three channels. The colorbar is shown on the right of the phase-directional image.

OCT devices and samples
A spectral-domain (SD-) OCT device with two slightly different configurations, denoted as low-resolution and high-resolution configurations, were used for the OCT volumes acquisition. The light source was a superluminescent diode (SLD, IW-M-D840-HP-I-CUS, Superlum, Ireland). The emitted light was split equally by a fiber coupler (TW850R5A2, Thorlabs Inc.) and incident into the reference and sample arms of a fiber Michelson interferometer. In the reference arm, a variable neutral density (ND) filter (NDC-50C-2M, Thorlabs) was used to control the reference power and a two-paddle polarization controller (FPC020, Thorlabs) was used to align the polarization state. In the sample arm, the light was collimated by a collimator (F280APC-850, Thorlabs) to a beam diameter of 4.0 mm, passed through a two-axis galvanometric scanner (GVS102, Thorlabs), and focused on a sample through an objective (LSM02-BB, Thorlabs, 18-mm focal length), so that the in-focus lateral resolution was 4.9 µm. The spectral interference signal was detected by a spectrometer which acquired the spectrum at a detection speed of 50,000 lines/s with 2048 pixels. The details of the spectrometer is described later. The spectrum was digitized by a Camera Link frame grabber (PCIe-1433, National Instruments, TX) and transferred to a PC.
For low-resolution configuration, the light source operated at a narrow-band mode with a center wavelength of 820 nm and a bandwidth of 70 nm. For this configuration, the spectrometer is a prototype device (Horiba, Kyoto, Japan) equipped with a line camera (spL4096-140km, Basler AG, Germany). The axial resolution and the depth-pixel separation were 4.2 µm and 3.1 µm in air, respectively. The sensitivity was measured to be 90 dB with a probe power of 2.8 mW on the sample.
For high-resolution configuration, the light source operated at a wide-band mode with a center wavelength of 840 nm and a bandwidth of 100 nm. The spectrometer was a Cobra OCT spectrometer (CS800-850/140-250-OC2K-CL, Wasatch Photonics, NC). The axial resolution and the depth-pixel separation were 3.1 µm and 2.5 µm in air, respectively. The sensitivity was measured to be 100 dB with a probe power of 6.0 mW on the sample.
A volumetric OCT was obtained by rescaling, numerical dispersion compensation, and fixedpattern-noise (FPN) removal. Here complex-median-line-subtraction algorithm [28] was used for FPN removal. Depth-independent bulk-phase error due to system imperfection and air turbulence was numerically corrected after reconstructing the OCT volume using the method presented by Oikawa et al. [29]. It is noteworthy that this bulk-phase-error correction is necessary for applying the numerical masks correctly. Figure 4 shows the spatial frequency spectra of an en face OCT image [ Fig. 4(a)] without [ Fig. 4(b)] and with [ Fig. 4(c)] the bulk-phase-error correction. Here the zero frequency is at the center. Without bulk-phase-error correction, the spatial frequency spectrum is significantly shifted along the slow scan frequency direction. Hence, the center of the spectrum and the center of the numerical masks described in Section 2 cannot collocate with each other. This shift is corrected by the bulk-phase-error correction, and hence, the numerical aperture mask can be properly applied. Four ex vivo bovine Achilles tendon and four chicken breast muscle tissues were used to validate the method. These samples were cut into a size of 10 mm × 10 mm. During the measurement, saline solution was put on the sample surface to avoid drying. The lateral measurement area was 1 mm × 1 mm and was sampled with 512 × 512 A-lines. Figure 5 summarizes the directional analysis of the bovine Achilles tendon. For this measurement, the high-resolution configuration of the SD-OCT was used. Figure 5(a) shows the angle distribution of the structural orientation of an en face OCT [ Fig. 5(b)]. The rotation angle is defined as the wheel diagram in which the fan-shaped opening is aligned vertically for 0 • , i.e., along the slow scan direction. The plot shows a distinct peak and the centroid of the peak was 60 • . So, the dominant fiber orientation is 150 • , i.e., 60 • + 90 • , which corresponds to the fiber orientation visible in the en face OCT image.    Figure 6 summarizes the directional analysis of the chicken breast muscle. For this measurement, the low-resolution configuration of SD-OCT was used. Intensity directional analysis was applied to the en face OCT [ Fig. 6(b)]. In the angle distribution of the structural orientation [ Fig. 6(a)], a distinctive peak (primary peak) is found at 159 • . In addition, a small but broad secondary peak is found at around 30 • to 90 • . In the phase-directional image perpendicular to the dominant orientation [ Fig. 6(d)], purple-green stripes are seen along the fiber structure. On the other hand, in the phase-directional image parallel to the dominant orientation [ Fig. 6(f)], fine purple-green stripes are appeared orthogonal to the fiber orientation. This may indicate the crimps of the muscle fiber (see Discussion in Section 5.2.2 for details). This muscle crimp would account for the secondary peak of the angle distribution in Fig. 6(a). The crimp pattern was found in all four chicken breast muscle samples we examined. On the other hand, all four bovine Achilles tendon samples including the one discussed in Section 4.1 did not show crimp patterns.

Possible projection artifact
In Section 4.2, we have interpreted the vertical stripe appearance in the cross-sectional phasedirectional images [Figs. 6(e) and (g)] as the vertically extending fiber structure. However, this appearance also can be interpreted as a projection artifact.
In our current interpretation model (Section 5.2), we assume that we are measuring the surface profile of a high refractive index structure in a low refractive index medium. On the other hand, the phase of OCT is sensitive to the path-length of the probe beam, which is proportional to the net refractive index in the path. And hence, the phase-directional signal in the deeper region can be affected by the refractive index of the superior tissue. It can cause a projection artifact, which is similar to that can be seen in Doppler OCT. To further validate this issue, histology investigation would be helpful in future.

Interpretation of phase-directional image by using numerical simulation
In this section, we present numerical simulations of the phase-directional imaging in order to interpret the results in Section 4. We first describe the system and sample models, and show the simulation results and discuss the interpretation of the experimental results.

System and sample models
For this simulation, the OCT was modeled as a simple Michelson-interferometer based SD-OCT. The light source has a Gaussian spectrum with a center wavelength of 840 nm and the bandwidth of 100 nm.
In this simulation, the sample is represented by a wavy surface structure z s (x, y) as depicted in Figs. 7(a) and (b). The surface can be considered as the interface between cylindrical fibers and surrounding medium. Here x and y are the lateral positions, z is the depth position and z = 0 represents the zero-delay depth. Note that we will add crimp structure to this cylindrical model later. In the particular geometry of Fig. 7(a), denoted as positive geometry, the sample surface is directed upward. In this geometry, z s (x, y)>0 and the probe beam illuminates the sample from the top. Although our experiments in Section 4 were performed with this positive geometry, we also performed the simulations with the "negative geometry" [Fig. 7(b)]. For this geometry, z s (x, y)<0 and the surface directs toward the zero-delay. The probe illuminates the sample from the top also in this geometry.
In order to imitate the fiber structures of the bovine Achilles tendon and the chicken breast muscle in more details, we added the crimp structure to the cylindrical fiber model as depicted in Fig. 7(c). And the sample surface z s (x, y) is written as where Cyl(x) is a repeating elliptical cylindrical surface along the y direction, where D x and D z are the diameters of elliptical cylinder along x and z, respectively. Each integer value of j = · · · , -2, -1, 0, 1, 2, · · · correspond to each cylinder. Here each cylinder represents each fiber in the sample. The sinusoidal term of Eq. (1) accounts for the crimping structure along the y direction, where A y and P y are the amplitude and the period of the crimp, respectively. T x and T y are bulk tilts of the sample along x and y directions. Z d is the constant depth offset to select positive and negative configurations. In this simulation, we set the size of the elliptical cylinder based on the measured size of the bovine Achilles tendon [ Fig. 7(c)] as D x = 105 µm and D z = 110 µm.
By assuming the probe beam is backscattered (backward-reflected) at the surface, the electric field of the probe beam at the spectrometer becomes while the reference beam isẼ where k is the wavenumber,S(k) is the spectrum of the light source, and z 0 is the common path-length of the probe and reference arms.
In the simulation, we numerically generated these two fields, summed them, and absolutesquared the sum to yield the spectral interference fringe. And OCT signal is generated by inversely Fourier transforming the fringe. The depth pixel separation, lateral scanning range and pixel separations are set to be the same with the high-resolution configuration (see Section 3). Phase-directional images are then created from this simulated OCT signal by the same method which has been applied to the experimental data (see Section 2.2).

Simulation results and interpretation
We first simulated the most simple example, a straight and flat fiber bundle with no crimp structure (A y = 0) and no bulk tilt (T x = T y = 0). For the simulation of positive geometry, Z d was set to be a large positive value to ensure z s (x, y)>0, while for negative geometry, Z d was set to be a large negative value to ensure z s (x, y)<0. The results are summarized in Fig. 8 (1), phase-directional images resembling the experimental results can be generated. To mimic the bovine Achilles tendon [ Fig. 5], we set A y = 0 µm, T x = 0 µm/mm, and T y = 15 µm/mm. This model represents a tilted straight fiber bundle. The resulting pseudo-color phase-directional volume is shown as rendered volumes in Figs. 9(a) and (b) with perpendicular and parallel sub-aperture sets, respectively. The perpendicular sub-aperture set results in purple-green stripes, while the parallel sub-aperture results in a homogeneously green appearance. Note that the sample model is tilted for y as indicated by yellow arrows in Fig. 9(b). Both of these simulation results well resemble the experimental results in Fig. 5. This resemblance indicates that the tendon sample would have a cylindrical fiber structure tilted along the fiber direction as the bottom-left of Fig. 5 is higher than the top-right.
To mimic the chicken breast muscle [ Fig. 6], we set A y = 20 µm, P y = 55 µm, T x = 0 µm/mm, and T y = 0 µm/mm. This model represents a crimping fiber bundle. The volume rendering results are shown in Figs. 9(c) and (d) with a perpendicular and parallel sub-aperture sets, respectively. The perpendicular sub-aperture sets generated straight purple-green stripes along the fiber direction, and it well resembles the experimental result [ Fig. 6(d)]. The parallel sub-aperture sets  results in high-frequency alternation of purple and green along the fiber, it also well resembles the experimental result [ Fig. 6(f)]. So, this result suggests that the chicken breast muscle would be a crimping fiber bundle.

Visibility of sub-OCT-resolution structures
In order to discuss the visibility of sub-OCT-resolution structures in the phase-directional image, we performed another set of simulations as shown in Fig. 10. Here the simulated structure is the same with Figs. 9(c) and (d) (chicken breast muscle with D z = 110 µm, P y = 55 µm) except for the amplitude of the crimp and the tilts (T x and T y ). For these simulations, the height of the crimp (A y ) was set to 1 µm, i.e., the peak-to-valley crimp height is 2 µm, which is smaller than the depth resolution of the OCT (3.1 µm). The tilts were set as T x = T y = 0 µm/mm (no tilt) for Figs. 10(a) and (b), T x = 1 µm/mm, T y = 0 µm/mm (horizontal tilt) for Figs. 10(c) and (d), T x = 0 µm/mm, T y = 1 µm/mm (vertical tilt) for Figs. 10(e) and (f), and T x = T y = 1 µm/mm (diagonal tilt) for Figs. 10(g) and (h).
In all configurations, the sub-OCT-resolution crimp structure is clearly visible as a horizontal red-green stripe pattern [Figs. 10(b), (d), (f), and (h)]. Although the tilts affect the image appearance, it does not have a significant impact on the visibility of the red-green stripes.

OCT intensity image and phase-directional image
Although the similar fiber structures are visible in both the intensity OCT image and the phasedirectional image, they are based on different properties of the sample. The intensity OCT shows the scattering property of the sample. On the other hand, the phase-directional image represents the intra-tissue surface structure (as discussed in Section 5.2), optical-path-length information, or similarly refractive index inhomogeneity in the sample.

Relationship between the first and second steps
In our previous work, we have used a complex OCT signal, i.e., not an intensity only or phase only signal but a fully complex signal, for directional processing [27]. On the other hand, in the present study, we used two-step evaluation, the first step uses only the intensity signal and the second step uses the phase information. This two-step approach was selected for the following reason.
The purpose of the entire analysis is to evaluate the lateral and depth oriented property of the tissue structures. The first step is to determine the en face orientation of the structure-of-interest in the tissue, which is visible in the OCT intensity image. So, it is reasonable to use the intensity OCT.
On the other hand, the second step is to determine the depth-oriented sub-OCT-resolution property of the structure-of-interest. So, it is mandatory to use the phase information. In addition, this depth-oriented property has not been encoded in the intensity but solely in the phase. And hence, we have used the phase-only signal, not the full complex signal.

Conclusion
We demonstrated computational multi-directional OCT, which visualized the directional property of tissue microstructure, and realized through post-measurement signal processing of a complex OCT volume. This signal processing consists of two steps. In the first step, the intensitydirectional analysis, the dominant fiber orientation in the en face plane was determined. In the second step, the phase-directional imaging, pseudo-directional image was generated. This processing method was applied to bovine Achilles tendon and chicken breast muscle samples, which showed clear en face stripes. By comparing with the numerical simulation results, the bovine and chicken samples were inferred to have a tilting straight fiber bundle structure and crimping fiber bundle structure, respectively. The experimental results and their simulation-based interpretation suggest that our method is useful for microstructural analysis of musculotendinous tissues.
Funding. Japan Science and Technology Agency (JPMJMI18G8); Japan Society for the Promotion of Science (18H01893, 21H01836).