2.1.

Tissue

Two euthanized adult rats were obtained from the tissue sharing program of Research Animal Resources at the University of Minnesota. Brains were dissected and fixed with 10% buffered formalin for 72 h before imaging.

2.2.

Serial Optical Coherence Scanner Imaging

The imaging system and associated experimental procedures have been described in detail by Wang et al.¹⁶ In brief, SOCS integrates an MC OCT and a microtome slicer to accomplish the reconstruction of macroscopic tissues. Built on a polarization-maintaining-fiber-based polarization-sensitive OCT in the spectral domain,²⁵ the MC-OCT characterizes anisotropic tissues by the optical property of birefringence and incorporates the Doppler flow information when applicable. The light source, a superluminescent diode operating at 840-nm central wavelength with a 50-nm bandwidth, yields an axial resolution of $5.5\mu \mathrm{m}$ in tissue (with a refractive index of $\sim 1.4$). The lateral resolution determined by a scan lens ($f=36\mathrm{mm}$) is about $15\mu \mathrm{m}$. Interference of the backscattered/reflected light from the sample and reference arms is detected by a home-built spectrometer that simultaneously acquires spectra on two orthogonal polarization states. Inverse Fourier transform of the spectral oscillations in $k$-space yields complex-valued depth profiles (A-line) for each polarization channel, which can be represented in polar form as $A_{\mathrm{1,2}}(z)\exp \{i{\varphi }_{\mathrm{1,2}}(z)\}$, where $A$ and $\phi$ denote the amplitude and phase, respectively, along depth $z$, and the subscripts 1 and 2 correspond to the cross-coupled and main polarization channels. The contrasts for ex vivo brain imaging were generated from the amplitude and phase of the complex-valued depth profiles. These include reflectivity $R(z)\propto A_1{(z)}^2+A_2{(z)}^2$, cross-polarization $C(z)\propto A_1{(z)}^2$, retardance $\delta (z)=\arctan [A_1(z)/A_2(z)]$, and optic axis orientation $\theta (z)=[{\phi }_1(z)-{\phi }_2(z)]/2$. The optic axis orientation, however, is a relative measurement, bearing a time-variant offset that needs to be removed. The offset originates from an arbitrary delay between the optical paths of the two PMF channels. To obtain the absolute optic axis orientation, a retarder film was included as a calibrating reference and imaged together with the tissue. Construction and operation details of the imaging system can be found in our previous publications.^16,25

The brain sample was mounted on a vibratome slicer positioned under the scanner. A volumetric scan (optical section) contains 300 cross-sectional frames (B-line) with 1000 A-lines in each frame, covering a field of view of $7\times 7\times 1.78{\mathrm{mm}}^3$ ($xyz$) with a corresponding voxel size of $7\times 23\times 3.47\mu {\mathrm{m}}^3$. After imaging one optical section, a superficial slice is removed allowing deeper regions to be imaged. The procedure is repeated until the whole sample is imaged. As the tissue is kept unmoved during the imaging procedures, 3-D reconstruction of the entire sample is achieved by stacking the slices without requiring complicated registration algorithms.¹⁶

Two rat brains were scanned, one sectioned in sagittal planes (shown in Figs. 2Fig. 3–4, 6, and 7) and the other in coronal planes (shown in Fig. 5). The sagittal sections consist of 28 slices of 200 μm each, and the coronal sections consist of 66 $100\text{-}\mu \mathrm{m}$ slices.

2.3.

Image Reconstruction

The 2-D en-face images unveil surface and subsurface features within a section. The volumetric datasets of reflectivity, retardance, and optic axis orientation contrasts were projected onto the $xy$-plane for each optical section. The pixel intensities in en-face reflectivity and retardance images were computed by taking the corresponding mean value along the depth ($z$) direction. The pixel intensities in en-face orientation images were determined by the peak of a histogram formed by binning the measured orientation values along the depth into 2 deg intervals. The depth range used in calculations matches the physical slice thickness. 3-D images of the entire sample block were reconstructed in two ways: by stacking the 2-D en-face images or by stitching the volumetric datasets of all sections. The en-face stack facilitates the global structure identification at a mesoscopic resolution ($15\times 15\times 100\mu {\mathrm{m}}^3$, the resolution in $z$-axis is constrained by the slice thickness). On the other hand, stitching the volumetric datasets preserves the microscopic resolution as recorded ($15\times 15\times 5.5\mu {\mathrm{m}}^3$). For optimal fusing outcome, the starting point along the depth was adjustable to match the ending point of the previous section. The cross-polarization images were used for the 3-D stitching, as the contrast describes the nerve fibers and preserves decent signal intensity in deeper regions. The details of the 3-D reconstructions can be found in our previous publication.¹⁶

The reflectivity image portrays the anatomical structures including the fiber tracts, and the retardance and cross-polarization images particularly probe the nerve fibers. They are used as an original dataset to apply ST. The optic axis orientation contrast provides quantification of the in-plane orientation of nerve fibers; therefore, it is used to compare with the results of the ST analysis. Details of the comparison methods are described in the following sections.

2.4.

Structure Tensor Analysis

The ST is a second-moment matrix, which is computed from the gradients of the image data.²⁶ It describes the dominant direction in a local neighborhood of a specific point. The computations are performed based on a pixelwise method.^20,24 The procedures are described as follows:

Noise, especially speckle, is a nontrivial problem in coherent and low-coherent imaging systems. Because the speckle pattern in SOCS images can be sensitively captured by the ST analysis, its removal by appropriate filtering is important. We used a nonlinear anisotropic diffusion filter described by Kroon and Slump.²⁷ The filter was applied as a preprocessing step on the original SOCS images. The effectiveness of filtering is evaluated in Sec. 3. The ST algorithms were implemented in MATLAB using customized codes.

2.5.

Comparison of Computed and Measured Orientations

The 2-D computational orientation maps are compared with the en-face optic axis orientation images obtained by SOCS. The orientation differences in selected ROIs are inspected. Since the inclination angle is not measured optically, the comparison of computational and measured orientations is not directly available for 3-D. Instead, we projected the computed 3-D vector onto the $xy$-plane and created an en-face image for each optical section using a histogram analysis, as described in Sec. 2.3. The 2-D images derived from the 3-D ST analysis are then compared with the en-face optic axis orientation images.

2.6.

Tractography

Since the ST matrix has the same construction as the tensor data in DTI, tractography tools used for the dMRI technique are readily applicable. We conducted tractography on the 3-D ST data using the Diffusion Toolkit.²⁸ As the tracts are created based on the eigenvectors with respect to the greatest eigenvalues, while in ST the fiber orientation is represented by the vector corresponding to the smallest eigenvalue, we replaced the eigenvalues (${\lambda }_i$) with $1/{\lambda }_i$ and rebuilt the ST matrix. We demonstrated the results with the interpolated streamline algorithm,²⁹ while application of the fiber assignment by continuous tracking method,³⁰ the second-order Runge Kutta²⁹ and the tensorline³¹ algorithms showed similar pathways. The tracts are visualized in TrackVis.²⁸ For detailed inspection, ROIs are selected and fibers passing through the nodes are examined. The pipeline of the entire processing is shown in Fig. 1.

Fig. 1

Processing pipeline of structure tensor (ST) analysis of serial optical coherence scanner (SOCS) images. The optical contrasts are generated for each volumetric scan. En-face images of reflectivity (R) and retardance ($\delta$) are calculated for two-dimensional (2-D) ST computation. Stacked en-face retardance images and stitched optical sections of the cross-polarization (C) contrast are utilized for three-dimensional (3-D) ST computation. The computed orientation maps are obtained through eigen-decomposition of ST and validated by the en-face optic axis orientation images. Tractography is conducted based on the 3-D ST.

3.1.

En-Face Representations of Structure Tensor

We applied ST on various contrasts of en-face images. A nonlinear diffusion filter was applied as a preprocessing step on the original images before the ST computation, and the pixel size was interpolated to be $7\times 7\mu {\mathrm{m}}^2$. The results with the reflectivity and retardance images are shown in Fig. 2. The left column shows the original images. The reflectivity image is the traditional OCT measure describing the gross structure, but the white matter can be either brighter or darker than the gray matter depending on the fiber orientation with respect to the illumination beam. The retardance image is generated from polarization-sensitive measurement and especially targets the white matter where the birefringence property of myelinated fibers alters the polarization state. The right column shows the computed fiber orientation maps. The images share the same color-coding given on the color wheel, while the brightness is determined by the original SOCS images (left). The results indicate that the computed orientations of the fiber tracts are consistent between the two contrasts. The best agreements are found in regions where fiber identifications are highlighted on both contrasts indicating in-plane alignment of those fiber bundles (e.g., the labeled fiber groups around the medial part of the thalamus). Some discrepancies are observed in the fiber clusters caudal to the thalamus, attributed by the different characteristics revealed on the original images. The result also suggests that the ST analysis can be applied to the images of conventional OCT where the polarization information is not available. As the retardance contrast provides the most robust identification of the white matter on the en-face images,^14,17 we apply the 2-D ST on this contrast for further evaluation and comparison.

Fig. 2

The 2-D ST analysis on (a) en-face reflectivity and (b) retardance images. Left: original images, right: orientation maps computed by eigenvector of ST. Color wheel shows the orientations. Brightness is controlled by the pixel intensities in corresponding en-face images. sm: stria medullaris thalamus, f: fornix, mt: mammillothalamic tract, and fr: fasciculus retroflexus. The round area in blue indicates the ROI 4 for Fig. 4(b).

ST is sensitive to abrupt changes in image intensity; as a result, noise reduction plays an important role in obtaining smooth and accurate representations of fiber orientations. Therefore, we applied a nonlinear diffusion filter before ST computation and compared the results with those without any filtering (Fig. 3). Two sizes of smooth kernels, $K_{\mathrm{10,2}}$ and $K_{\mathrm{20,4}}$, were used for the tensor matrices. $K$ denotes a Gaussian kernel with the first subscript being the number of neighboring pixels ($7\mu \mathrm{m}/\text{pixels}$) included and the second subscript being the standard deviation of the Gaussian kernel. The maps are accompanied by the histograms of fiber orientations in an ROI (dashed box). The ROI is placed in the optic tract where coherent orientation is expected in the fiber bundles. When no filter was utilized on the original image, the orientation map is hardly immune from noise (top left) for a smoothing size of $28\mu \mathrm{m}$ ($K_{\mathrm{10,2}}$), and the histogram presents a scattered orientation distribution in the optic tract. Increasing the size of the smooth kernel and the weight of the neighboring points enhances the consistency of orientation representation (right column), but results in degraded resolution and misidentification of small crossing fibers if existing. To minimize the problem, a preprocessing step with a nonlinear diffusion filter dramatically reduces the noise while keeping the spatial resolution less affected. We empirically determined to use a smoothing kernel size of 30 to $50\mu m$ with the preprocessing filter in the following sections.

Fig. 3

Filtering and kernel size effects on the computational orientation maps: (a) without filter; and (b) with nonlinear diffusion filter applied on SOCS images. The sizes of the smoothing kernels are $K_{\mathrm{10,2}}$ (left) and $K_{\mathrm{20,4}}$ (right), respectively, where the first subscript is the number of neighboring pixels included in smoothing and the second subscript is the standard deviation of the Gaussian kernel. The color-coding is in Fig. 2. The white boxes indicate the region of interest (ROI), whose fiber orientation distributions are plotted with histogram (white bars) on individual images (range: $\pm 90\mathrm{Deg}$, binning: 2 deg).

3.2.

Comparison with Optical Measurement

To evaluate the 2-D ST performance, we compared the computational orientations with the optic axis orientation measurements. Figure 4(a) shows a sagittal slice of the computed and measured orientation maps. The orientation values share the same color coding given by the color wheel. The brightness of the pixels is controlled by the en-face retardance values. The orientation maps demonstrate a remarkable agreement with less than a 10 deg difference in well delineated fibers where higher brightness (retardance) is usually seen. These areas include the internal capsule (ic), the stria medullaris thalamus (sm), the cingulum (cg), and the fiber tracts in the caudate putamen (CPu). More deviations are observed in the white matter regions with low brightness.

Fig. 4

Comparison of computed and measured orientations: (a) The orientation maps of a sagittal section share the same color-coding (color wheel), and the brightness is controlled by the en-face retardance values. CPu: caudate putamen, TH: thalamus, cc: corpus callosum, and ic: internal capsule. (b) Mean (circles) and standard deviation (bars) of the difference between computed and measured fiber orientations are plotted for all ROIs and respective slices. Pixels having low retardance ($<233\mathrm{nm}$) are excluded in the statistical analysis. ROIs of the first three plots are indicated on (a). ROI of the last plot contains multiple fiber bundles as shown in Fig. 2 including f, fr, mt, and sm. The ROIs are across multiple slices [$x$-axis in (b)].

We then performed quantitative comparisons in specific areas of the white matter [Fig. 4(b)]. The ROIs 1, 2, 3 [indicated by the boxes on Fig. 4(a)] contain small fiber tracts in the CPu, dense fiber tracts in the ic, and fiber tracts in the lateral thalamus, respectively. The ROI-4 contains fiber groups around the medial thalamic region including the stria medullaris thalamus (sm), the fornix (f), the mammillothalamic tract (mt), and the fasciculus retroflexus (fr), as indicated in Fig. 2. The ROIs are across multiple slices where the selected architectures are clearly visible. Figure 4(b) shows the plots of mean and standard deviation of the difference between the computed and measured fiber orientations. Pixels having retardance values of 233 nm and higher are included in the statistical analysis. The orientation differences ($\text{mean}\pm \text{standard deviation}$) are $-0.2\pm 15.5\mathrm{deg}$, $-5.7\pm 22.1\mathrm{deg}$, $-14.6\pm 28.7\mathrm{deg}$, and $-7.9\pm 29.41\mathrm{deg}$, respectively, for ROIs 1, 2, 3, and 4. The minimal differences indicating good agreement are seen where individual fiber tracts are clearly identified (ROI-1). In addition, the computational orientations are viable in fiber bundles where the dominant features can be captured (ROIs 2 and 4). The orientations are less comparable with greater variations observed in the lateral thalamus (ROI 3) where fiber identification becomes less reliable.

3.3.

Three-Dimensional Orientation Maps

The 3-D ST achieves quantification of the fiber orientation in the whole brain space, which fills the gap of current 2-D orientation measurements by SOCS. First, we conducted the ST analysis on the en-face stack of the retardance images. The voxel size of the data was interpolated to be $25\mu \mathrm{m}$ isotropic. This dataset provides a global view of the fiber organization and captures the neural fibers running through the $xy$-plane with an inclination angle. Figure 5(a) demonstrates the orthogonal views of the 3-D orientation map derived from the ST. The fiber directions are indicated by the color-coding on each map, and the brightness is controlled by the retardance value. The coronal view (left) shows various fiber groups including the anterior commissure (ac) connecting the two hemispheres in the left-right direction, the f aligned in the superior-inferior direction, and the fibers in the CPu running in the rostral-caudal direction. Color transitions are seen in the fimbria of the hippocampus (fi) and the posterior branch of the anterior commissure (acp). Similar results are seen on the sagittal (middle) and horizontal (right) sections as well. The orientations are consistent with the prior knowledge from dMRI results.^32,33 Figure 5(b) demonstrates a microscopic diffusion tensor atlas (resolution: $50\mu \mathrm{m}$) of the Wistar rat brain.³² Orthogonal views of 3-D orientation maps are shown at similar locations. The color-coding for orientation is the same as Fig. 5(a). The majority of the fiber orientations obtained by ST are well correlated with the diffusion tensor images. The major deviations in the ST results are observed in the corpus callosum (cc) and external capsule (ec) which form large white matter areas in the absence of the detailed features of fiber bundles. The lack of a dominant direction in local features may induce errors in the orientation estimations within large white matter regions.

Fig. 5

(a) The 3-D ST for fiber orientations at mesoscopic resolution. The ST is calculated from a stack of en-face retardance images and the orthogonal images are visualized in Vaa3D³⁴ (left: coronal, middle: saggital, and right: horizontal). The color represents the fiber orientations (red: left-right, green: anterior-posterior, and blue: superiorinferior). Directions are labeled on the individual viewing plane. The image intensity is masked by the retardance value which highlights the white matter. CPu: caudate putamen, acp: posterior branch of the anterior commissure, cc: corpus callosum, ec: external capsule, f: fornix, and fi: fimbria of the hippocampus. Sample size: $7\times 7\times 5.5{\mathrm{mm}}^3$. (b) Microscopic diffusion tensor atlas of a Wistar rat brain (resolution: 50 μm) shows the correspondence with the same color-coding. The fiber orientation maps are shown superimposed on the fractional anisotropy images, which were obtained from the public online dataset at the Duke Center for In Vivo Microscopy.

The 3-D dataset generated by stitching the volumetric optical sections of serial scans provides higher resolution for delineating the fiber map. ST analysis applied on this dataset provides a favorable solution for both large-scale orientation quantification and comprehensive connectivity investigations. To achieve more efficient computation, the dataset was downsampled to a voxel size of $25\mu \mathrm{m}$. Figure 6 shows a 3-D orientation map of the right hemisphere of a rat brain. The colors represent the orientations as follows: red: rostral-caudal, green: left-right, and blue: superior-inferior. The brightness is controlled by the pixel values of the cross-polarization images. Figure 6(a) illustrates the orthogonal views demonstrating fiber tracts with various directions and patterns. Directional changes of fiber tracts can be sensitively detected. On the coronal plane, the neural tracts through the putamen are aligned more horizontally in the lateral region (ROI 1) and vertically in the upper region (ROI 2). Similar color transitions are shown on the horizontal and sagittal planes. Moreover, results obtained from the 3-D ST analysis clarify the ambiguity of fiber alignments on 2-D ST images. For example, the orientation map on the horizontal plane shows that the fiber bundles in the ic run along the superior-inferior direction rather than the left-right direction which would be the apparent intuition from the 2-D plane. Fiber bundles aligned in the left-right direction are less visible, probably due to the fact that the optical systems barely detect the fibers along the illumination beam. This yields low backreflection and decreased signal intensity compared with the surrounding gray matter; therefore, some of the fiber bundles are not captured. Volume rendering of the 3-D orientation map is implemented in Vaa3D³¹ and shown in Fig. 6(b).

Fig. 6

The 3-D ST analysis on high-resolution reconstruction of rat brain. The dataset is generated by stitching the optical sections of the cross-polarization contrast. (a) Orthogonal views of the 3-D fiber orientations. ic: internal capsule. ROIs 1 and 2 indicate the fibers oriented along horizontal and vertical directions, respectively. The color sphere codes the orientations as red: anterior-posterior, green: left-right, and blue: superior-inferior. Directions are labeled on the orthogonal planes as well. Brightness is controlled by the dataset (cross-polarization). (b) Volume rending of the 3-D orientation map provides a perspective view.

The 3-D orientation obtained by ST analysis is not directly comparable with the optical measures as the optical axis orientation contrast in SOCS is restricted in 2-D. Therefore, we projected the 3-D computational orientation vectors onto the $xy$-plane and then generated en-face images to compare with the en-face optic axis orientation images. Pixel values on the en-face computational orientation image are determined by a histogram approach similar to that used for en-face optic axis orientation. Figure 7 shows the en-face images of a sagittal section. The computation (left) and the measurement (middle) are well correlated in most of the white matter, as shown by the absolute difference between the two images (right). Fiber orientations obtained by the two approaches closely match with the differences less than 10 deg in most regions where the image intensity is high and neural fibers are clearly traceable. These regions include the ic, the tracts in the putamen, the optic tract, the superior thalamic radiation, and the local fibers within the thamalus. More deviations are seen in the fiber bundles going through the plane.

Fig. 7

The 3-D orientation by ST analysis is (a) projected onto the $xy$-plane and (b) correlated with the en-face optic axis orientation. The color wheel indicates the orientations. The brightness is controlled by the en-face retardance values. The absolute difference of computed and measured orientations for pixels with retardance greater than 233 nm is plotted on the en-face retardance image (c). The color bar represents a narrow range (0 deg to 35 deg) for visualizing the mismatch.

3.4.

Tractography

Tractography can be readily performed using the tools designed for diffusion MRI data. Figure 8 shows a representative tractography based on the ST analysis applied on the en-face retardance stack ($7\times 7\times 5.5{\mathrm{mm}}^3$). The tracts are color coded by their local orientations according to the directions shown in the cube (bottom-right corner). The original dataset is overlaid in grayscale for a better interpretation of the anatomy. Figure 8(a) illustrates the pathways beneath the cc. Only 14% of the tracts are shown due to the high density. Tracts in the mass white matter areas such as cc and ec are excluded, because the orientation estimation of the 3-D ST on these shell-like geometries does not align with the direction of its constituting fibers. To explore the fibers passing through a specific region, a sphere ROI indicated by a black arrow on Fig. 8(b) is placed at the conjunction of f and ac. The localization of the ROI is facilitated with the visualization of 3-D anatomical images of the en-face retardance stack. The tracts on Fig. 8(b) indicate that the fibers passing through this region primarily include f, fi, and ac. The ac is further branched in the anterior and posterior directions (aca and acp). Tracking of the ac is better visualized when the ROI is moved more laterally [Fig. 8(c)].

Fig. 8

Tractography based on SOCS images. Tracts in color are computed from ST applied on the en-face stack of retardance. Colors in the cube represent the directions. Tracts are overlaid on the original data (grayscale). (a) Tracts excluding the cc and ec. Considering the high density, only 14% of the tracts are shown. (b) Fibers passing through a ROI at the junction between the ac and the f (gray sphere indicated by black arrow). (c) Fibers passing through another ROI (blue sphere) on the ac.