1. Introduction

Collagen is an essential protein in the human body that provides both structural and functional support in human organ tissue. Collagen forms as hierarchical bundles of fibers that begin at a molecular scale and end at the scale of entire organs [1]. On the macro-scale (millimeter to centimeter lengths), collagen fiber bundles coalesce and organize according to tissue function. This function can be to facilitate mechanical motion like in the case of tendon [2,3], provide a path of propagation for electrical signals in the heart or uterus [4,5], or simply support the curvature of the cornea [6,7]. To better understand and model these complex functions requires correspondingly detailed and complex models of the collagen architecture. Such models can be challenging to obtain in cases when collagen fibers form into complex three-dimensional networks, are embedded in muscle cells and tissue heterogeneities, and exhibit wide intra-specimen variance.

Optical Coherence Tomography (OCT) is a volumetric imaging modality with micron resolution capable of investigating tissue structures like collagen fibers in human and animal tissue [8]. Compared to similar optical microscopy techniques, OCT offers larger fields-of-view (0.5 cm) and can be used in-vivo. The combination of high resolution and field-of-view makes OCT an ideal candidate for studying collagen macro-structure as previously described in the anterior cruciate ligament (ACL) [9], cervix [10], and myocardium [11,12]. Modeling collagen fiber networks with OCT over many centimeters of tissue requires automated analysis of massive datasets. Image processing methods have been proposed that automatically extract fiber orientation from en-face OCT images using techniques based on image gradients [10,13] and the Radon Transform [14]. Numerous methods have been demonstrated for other imaging modalities like Second Harmonic Generation Microscopy (SHG) and Fluorescence Microscopy [15–20]. These microscopy techniques provide high contrast images of the microscopic organization of collagen and the extra-cellular matrix through fluorescent labeling and biomarkers that associated image processing algorithms inherently leverage in detection and modeling tasks. As a result, these processing methods tend to translate poorly to OCT images with lower signal-to-noise and spatial resolution.

While 2-D analysis has revealed valuable insight into the collagen structure-function relationship, in some cases 3-D analysis is necessary. Polarization-Sensitive OCT (PS-OCT) extracts additional information which can be used to quantify the orientation of collagen fibers. PS-OCT has been used to create 3-D tractography of collagen fibers in mouse hearts [21], the oral mucosa [22], cartilage [23], and neural pathways of the rat brain [24], but has seen limited demonstrations outside these models. Gan, et al. was the first to show that using structure information from Spectral Domain OCT (SD-OCT) images alone was sufficient to model 3-D fiber networks in the case of the human myocardium [12]. A limitation of this approach, however, is that it relies on the assumption that collagen fibers are roughly parallel to the tissue surface and is therefore incapable of modeling fiber networks with higher 3-D complexity. Furthermore, 3-D methods that work well for microscopy images like second harmonic generation microscopy tend to extend poorly to OCT images due to the unique challenges of image quality (high noise, low contrast, artifacts, etc$\ldots$). New approaches in deep learning have also been introduced for extracting fiber networks from fluorescence microscopy image volumes [25]. These methods produce promising results at high deployment speeds, but require massive amounts of data and experts for labeling training data.

More sophisticated 3-D methods are necessary to study particularly complex tissues like the human uterus. The human uterine wall comprises three distinct layers: the perimetrium, the outermost layer made of connective tissue; the myometrium, a thick middle layer consisting of smooth muscle cells and connective tissue including collagen and elastin [26]; and the endometrium, the innermost layer consisting of epithelial cells and stromal tissue. The smooth muscle cells comprising the myometrium are long and tubular; thousands of these cells align parallel to each other in sequence, forming a myocyte. Multiple myocytes further combine to form fascicles that are separated from each other by connective tissue composed mostly of collagen. This connective tissue acts as a sheath and provides structural scaffolding to which the muscle fibers can adhere. Our preliminary histology results have shown that in a pregnant uterus myometrium, smooth muscle cells and collagen account for approximately 80 and 20 percent of the total tissue volume, respectively.

In the uterus, collagen has a complex three-dimensional architecture that is associated with the remodeling process which occurs as a result of pregnancy. Previous studies using both DT-MRI [27] and histological processing have revealed organizational structure on the centimeter scale. Lutton, et al. were the first to do in-depth 3-D collagen fiber modeling of a human uterus specimen [28]. They used image processing of histology slides to create a 3-D model from a single tissue wedge. This approach revealed complicated organization at the micron scale, however, their method is too resource intensive to feasibly use in a multi-specimen studies. Therefore, a method capable of 3-D collagen fiber modeling that can both capture micron-scale structure and be used efficiently on multiple specimens is essential to improving our understanding of uterine tissue mechanics. Furthermore, little is know about typical uterine architecture so this method must also model fibers at random orientations with no prior assumptions of the tissue’s organization.

We propose a novel process for high speed imaging and analysis of 3-D collagen fiber organization using SD-OCT and an ensemble of image processing methods. Our approach has three phases: image enhancement, fiber orientation mapping, and fiber tractography. In this manuscript, we will detail the mechanics of each of these steps and validate them on simulated image volumes of fiber structures. The remainder of the manuscript will demonstrate the performance of our method on mosaic OCT image volumes of multiple human uterine tissue specimens. We show that the 3-D uterine architecture can be modeled in both pregnant and non-pregnant specimens and that quantitative analysis reveals new and exciting physiological structures and patterns such as fiber groupings.

2. Methods

A multi-phase image processing pipeline was developed to obtain quantitative fiber orientation information from stitched OCT volumes of both pregnant and non-pregnant uterus specimens. Figure 1 graphically shows this workflow to demonstrate the processing flow from image acquisition to 3-D fiber tractography. This section begins by explaining the 3-D geometry used to describe collagen fiber orientation and then describes the algorithms and implementations used in each phase: image enhancement, 3-D fiber orientation mapping, and 3-D tractography.

 

Fig. 1. Graphical representation of the 3-D fiber analysis pipeline. (a) Sample preparation, OCT volumetric imaging, and stitching to acquire mosaic volumes. (b) Image enhancement pipeline for reducing noise and artifacts while improving fiber contrast. (c) 3-D fiber orientation maps produced through a combination of b-scan and en-face analysis to obtain $\theta$ and $\phi$. (d) Histogram analysis of $\theta$ and $\phi$ reveals fiber groups while a 3-D particle filter technique is used to create fiber tracts.

Download Full Size | PDF

2.1 Fiber geometry

To obtain the 3-D orientation of collagen fibers from OCT image volumes requires two angles, $\theta$ and $\phi$. Figure 2(a) shows these how these angles are defined within the image volume coordinate system. The solid red arrow represents a single fiber of unit length which makes an angle with the positive Z-axis ($\phi$) and the positive X-axis ($\theta$). The angle $\theta$ can be directly visualized and measured from an en-face image; however, an orthogonal image plane likely does not exist to directly measure $\phi$. Instead, we measure its projections onto either transverse plane which we will refer to as $\phi _x$ and $\phi _y$ and are drawn as the dashed red lines in Fig. 2(a). The projections are simply the XZ or YZ image planes that intersect each collagen fiber. The fiber incline angle $\phi$ is calculated using the following trigonometric relationship:

 

Fig. 2. Coordinate system and fiber orientation diagram. (a) 3-D image volumes are imaged as a series of b-scans (XZ) with uniform spacing in Y to capture the entire uterine tissue slice. The solid red arrow represents a collagen fiber of unit length which is defined by the angle it makes with the positive Z-axis ($\phi$) and the positive X-axis ($\theta$). $\theta$ is measured directly from en-face image slices, while $\phi$ is calculated from $\theta$ and its projection onto either either transverse plane, i.e. $\phi _x$ or $\phi _y$ (Eq. (1)). (b) The measured $\theta$ and $\phi$ orientations are displayed in this manuscript according to the shown HSV colorwheels. Based on the orientation of the tissue during imaging, an angle of $\phi = 0^{\circ }$ (red) indicates fiber alignment perpendicular to the tissue surface while $\phi = 90^{\circ }$ (cyan) indicates fiber alignment parallel to the tissue surface.

Download Full Size | PDF

The color maps in Fig. 2(b) show how the angles $\theta$ and $\phi$ are displayed in this manuscript. The HSV colormap was chosen since it is cyclical, i.e. uses the same color to represent the minimum and maximum values (in this case $0^{\circ }$ and $180^{\circ }$).

With respect to the uterine tissue slice, $\theta$ can be interpreted as measuring a collagen fiber’s alignment within the tissue surface plane. $\phi$ can be interpreted as the fiber’s incline within the tissue slice. Because $\phi$ is defined with respect to the positive Z-axis, a measurement of $\phi = 0^{\circ }$ would indicate a fiber’s aligned vertically while a measurement of $\phi = 90^{\circ }$ would indicate a fiber’s aligned horizontally or parallel to the tissue surface. Defined anatomically, the Z-axis corresponds to the through-thickness dimension of the uterus sample while the X and Y -axes align with either the longitudinal or meridianal anatomical directions depending on location.

2.2 Image enhancement

Detection and analysis are image processing tasks that can be uniquely challenging to apply to OCT b-scans, particularly for the task of fiber orientation analysis. A custom image enhancement pipeline was developed (see Fig. 1(b)) with the specific goal of improving b-scan image quality in a way that produces the most accurate results from downstream fiber orientation analysis. Figure 4 outlines each step of this process.

To enhance the imaged collage fibers, we first confirmed the OCT features that correspond to collagen in uterine tissue. Generally in OCT images, collagen fibers are bright when compared to the signal background and surrounding tissue. Biologically, collagen is a fibrous structure that weaves between smooth muscle cell bundles. In an OCT image, this appears as a bright fibrous network surrounding lower intensity regions of muscle fibers. This structure is shown in Fig. 3, which provides a visual comparison between histology (Fig. 3(a)) and an en-face OCT image (Fig. 3(b)).

 

Fig. 3. Matching OCT and histology to verify tissue composition in a non-pregnant uterine sample. (a) Masson’s trichrome stain histology section from anterior location. The orange stain indicates smooth muscle cells (SMCs) and the blue stain indicates collagen fibers. (b) OCT en-face image from the matching uterine tissue sample. Visual comparison of (a) and (b) indicates that the bright intensity fibers in the OCT image are collagen and the darker intensity areas between the fibers are SMCs.

Download Full Size | PDF

The first objective for pre-processing was to improve image quality by de-noising and enhancing the contrast of fibrous tissue features. Speckle noise is the primary noise component in OCT images and a number of sophisticated techniques exist for de-noising OCT images while also preserving features of interest. For this application VBM3D was found to be most effective [29]. When applied to an image volume, VBM3D (a video-based variant of BM3D) uses a collaborative filtering approach to both de-noise and preserve smooth, connected features that span adjacent b-scans.

A troublesome artifact for fiber orientation detection is surface reflections. This artifact is common when imaging fresh tissue samples which have to remain hydrated during imaging and tend to have a thin layer of moisture on the tissue surface. When the surface of the tissue is flat and aligned perpendicular to the laser source, this surface liquid acts like a mirror which results in saturation of the a-line at that location. These bright vertical streaks can corrupt large portions of a b-scan and obscure the underlying tissue features. Both Fourier and image gradient based fiber orientation methods have no simple way to distinguish these artifacts from a vertically orientated collagen fiber so the artifact must be removed prior to analysis. We accomplish this using a simple wedge filter centered at the vertical orientation. While a vertically aligned fiber has some width and natural variation, reflection artifacts are perfectly vertical and have constant image intensity. This feature is leveraged by using a narrow wedge filter which will only filter lines that are exactly vertical without degrading image quality.

The final challenge is to improve image contrast. This challenge has two components which are the inherently narrow dynamic range of OCT images and the mean intensity fluctuations that occur from image stitching and depth-based intensity roll-off. This was corrected by first applying a homomorphic filter followed by histogram equalization. The homomorphic filter or "high-frequency emphasis" filter improves the sharpness of tissue and fiber edges while also creating an even mean intensity across the images. Intensity re-mapping is employed following the homomorphic filter to optimize the dynamic range. Pixel intensities on the range of $[0.138,0.333]$ were mapped to a normalized intensity range of $[0,1]$.

Figure 4 outlines each step of the pre-processing image enhancement pipeline an unprocessed b-scan from a non-pregnant (NP) anterior uterine tissue sample. In this example, reflection artifacts are prominent and obscure the underlying fiber architecture. Low contrast and speckle noise further obfuscate the tissue features. Figure 4(b) shows the same b-scan following artifact removal. Both the reflection artifacts and natural streaking texture of the b-scan have been smoothed out. Figure 4(c) is the denoised image acquired by applying VBM3D denoising to the image in Fig. 4(b). This step further cleans the image and creates smoother, more connected fiber tracts. Figure 4(d) is the final enhanced b-scan which results from homomorphic filtering of the denoised image in Fig. 4(c).

 

Fig. 4. Step-by-step process of b-scan image enhancement. (a) Un-processed b-scan from a non-pregnant (NP) anterior uterine tissue sample. (b) Resulting image from vertical reflection artifact removal. (c) VBM3D denoised image which smooths the image and fiber tracks by mitigating speckle noise. (d) The fully enhanced b-scan which results from homomorphic filtering of the denoised image. Image enhancement improves the Radon-based fiber orientation analysis. (e) Un-processed image patch from marked in (a) by the green square. (f) Processing image patch from the same location as (e). (g) Radon-based fiber orientation distribution for image patch (e). (h) Radon-based fiber orientation distribution for image patch (f). Reflection artifact removal has reduced the peak at $90^\circ$ while denoising and contrast enhancement have exaggerated the fiber-based peaks at $10^\circ$ and $60^\circ$. Scale Bar = 0.5 mm

Download Full Size | PDF

The image patches shown in Fig. 4(e) and Fig. 4(f) provide a closer look of the image enhancement effect. Both patches are from the same section of tissue marked by the red squares in Fig. 4(a) and Fig. 4(d). Figures 4(g) and 4(h) show the fiber distribution obtained by the Radon-based method when applied to the image patches in Fig. 4(e) and Fig. 4(f), respectively. The Radon-method (described in the next section) provides a distribution of the fiber orientation content of an image, such that the distribution peak corresponds to the dominant orientation of fibers in the image. The pre-processing has two appreciable effects on the fiber distribution. The first is that the reflection artifact filter significantly reduced the peak centered at $90^\circ$. This is the largest peak in the un-processed distribution (Fig. 4(e)) meaning that our automated method for determining fiber incline would falsely classify fibers in this patch as having a $90^\circ$ incline. The second effect is that the peaks corresponding to the true fiber orientation are enhanced. Comparing Fig. 4(e) and Fig. 4(f) which are displayed at the same scale, the peaks at $10^\circ$ and $60^\circ$ are taller and more distinguishable for the processed patch distribution which further improves the performance of the automated fiber incline detection.

2.3 B-scan fiber orientation mapping

Following pre-processing of the OCT image volume, the fiber incline orientation is iteratively measured in each b-scan and en-face orientation is measured in the perpendicular image plane (Fig. 1(c)). Fiber orientation analysis of stitched volumes requires segmentation of the tissue because only the orientation of the fibers should be analyzed and not the tissue shape. This requires edge detection of the top and bottom tissue surface so that the 3-D surface is acquired. First, mean a-line intensity is used to crop each b-scan laterally. The top surface is then detected for each a-line by extracting the max intensity location in the standard deviation filtered image. The bottom surface is detected by calculating the first pixel location in each a-line that hits the noise floor based on our imaging system’s typical intensity roll-off. The lower edge is additionally corrected to obtain a boundary that is smooth and varies with the upper edge. Applying this procedure to every b-scan obtains the 3-D tissue surface.

Despite fiber enhancement processing, automatically determining fiber orientation from b-scans using gradient based methods is infeasible due to the varying fiber width and the subtle vertical streaking that results naturally from normal OCT acquisition. To obtain accurate fiber orientation measurements a Radon-based method was employed that has been previously demonstrated on SHG and OCT images for ultra-fast fiber orientation analysis [14,30].

The Radon-based method uses a simple multi-step approach to obtain quantitative distributions of collagen fiber orientation within an image patch. First, the image patch is converted to Fourier space, bandpass filtered to pass fiber lengths of interest, and then shifted using fftshift() in MATLAB so the lowest frequencies are at the center of the image and frequency increases as you move radially away from the center. The qualitative orientation information captured in Fourier space is converted to a quantitative result by taking the Radon transform of Fourier space image. Due to the FFT shift, the center row of the Radon transform then corresponds to the concentration of fibers oriented at a given angle. The peak of this row is analyzed to extract the dominant fiber orientation within the image patch. Through manual measurement of the uterine b-scans, a conservative estimate of fiber diameter and length was determined and the band-pass filter range was accordingly set to [0.255, 4.59] $mm^{-1}$.

The free space setup of the OCT system used in this study resulted in b-scans which capture both the tissue sample and its surface topology. Consequently, we found drastic changes in the tissue surface topology to be commonplace which further complicated the task of using b-scans to sample and analyze the uterine fiber architecture. We addressed this challenge by using a custom segmentation procedure that we will refer to as "patch processing" where local patches of the tissue are segmented and analyzed iteratively rather than analyzing all the tissue at once. The "patches" are square regions that are generated parallel to the local tissue surface. An optimization procedure is used to find the patch size which fills the most space between the top and bottom tissue surface. Given a minimum patch size and fixed overlap between adjacent patches, the uterine tissue can be segmented and analyzed for arbitrary tissue topology. Patches are generated for each b-scan of the mosaic volume, the Radon method is applied to each patch, and the dominant fiber orientation within the patch is recorded as the distribution peak. This process is repeated for all b-scans in the image volume to produce a continuous fiber orientation map $\phi _{x,y}$ using 2D nearest neighbor interpolation.

Figure 5 shows the b-scan fiber analysis process for a single b-scan (Fig. 5(a)) from a PG sample as an example. The uneven surface of the tissue required a piece-wise segmentation procedure to measure $\phi _x$. This procedure involves first using edge detection to obtain the upper and lower tissue surfaces, the results of which are highlighted in red in Fig. 5(a). Next, the b-scan is sampled into overlapping square "patches" which are generated parallel to the local tissue surface. Each of these patches (shown as blue squares in Fig. 5(a)) are processed using the Radon-based fiber orientation method to obtain the dominant orientation within the patch.

 

Fig. 5. Example of the b-scan pre-processing pipeline with edge detection and image patching. (a) Enhanced b-scan from pregnant (PG) fundus specimen (Scale bar = 0.5 mm). Upper and lower tissue edges are drawn in red and processing patches are displayed as the blue squares. (b) Interpolated fiber orientation across the full b-scan. Red circles indicate the fiber orientation measured in each blue patch. (c) Inset of the green patch in (a) (Scale bar = 200 $\mu$m). (d) Fiber orientation distribution for the image patch in (c). The peak at $20^{\circ }$ corresponds to the dominant fiber orientation captured in the image patch.

Download Full Size | PDF

A magnified view of the green patch in Fig. 5(a) and its corresponding fiber orientation distribution is shown in Fig. 5(c) and Fig. 5(d). The distribution has a peak at $25^{\circ }$ (red asterisk) which matches the dominant orientation of fibers in that patch. Figure 5(b) shows the resulting fiber orientation incline graph for the full b-scan. The red circles are the dominant angles measured in each patch and the blue line is the interpolated result. By repeating this procedure for each b-scan in an image volume, we obtain a full $\phi _x$ fiber map which gives the projected fiber incline at each $(x,y)$ location of the imaged tissue.

The patches are constrained to be square for several reasons. One benefit of square images is that no corrections are needed to adjust for the difference in dimension size which is a necessary step for applying the Radon-based fiber distribution analysis to non-square images. This reduces the complexity of the algorithm and improves the computational speed. Rectangular patches could improve the sampling results by allowing the patch length to freely adjust to the exact tissue surface. It is non-trivial to adaptively adjust the patch width in this way and it’s unclear the extent to which this change would improve the accuracy of the fiber orientation analysis, therefore, the straightforward square patch approach was chosen.

OCT b-scans were acquired with a 15 $\mu$m lateral resolution and 6.5 $\mu$m axial resolution. This results in a non-square pixel size which means angle measurements will be distorted if a correction is not applied. Therefore, the $\phi _{x,y}$ maps are adjusted following interpolation according to the following transformation

2.4 3-D fiber orientation mapping

A 3-D description of a collagen fiber’s orientation is achieved by obtaining the associated angles $\theta$ and $\phi$ (Fig. 1(c)). The en-face fiber orientation $\theta$ is calculated using a traditional gradient-based method which has been previously demonstrated on OCT images of human heart tissue [13]. The gradient method takes an en-face OCT image as input, divides the image into small windows, and then calculates the fiber orientation detected within the windowed area. Local orientation is determined by the direction along which the image gradient within each window is maximized. The image gradient is obtained using 2-D, 3 x 3 pixel, horizontal and vertical Sobel filters $G_x$ and $G_y$, and the gradient direction is $\theta '=arctan(G_y / G_x)$. A D’Agostino-Pearson $\kappa ^2$ (normality) test was used to determine if a windowed image area contained a valid fiber orientation. A $\kappa ^2$ threshold value of 0.02 and window size of 51 x 51 pixels was selected for all analysis. Variance in tissue topology also required partial segmentation of the en-face images to calculate $\theta$. This was done by using the tissue surface obtained by edge detection during the b-scan fiber analysis and finding the pixels between the top and bottom tissue surfaces for each en-face plane. If any portion of a window contained pixels beyond the extent of the volume, than that window’s result was omitted.

An en-face orientation $\theta$ was calculated at every valid location of the OCT volume by iteratively applying the described gradient method to each en-face OCT image. The result is a 3-D map $\theta (x,y,z)$. The final step of the mapping process is to obtain the corresponding 3-D map of fiber incline orientation $\phi (x,y,z)$. This is calculated using Eq. (1), where $\phi _{x,y}(x,y)$ is assumed to be constant with depth.

2.5 3-D tractography

In addition to creating fiber orientation maps, we are interested in modeling the collagen architecture in 3-D to create a visual representation of the tissue (Fig. 1(d)). To this end, we employed a tacking technique called particle filtering [31]. The goal of particle filtering in this application is to find the most likely path through the OCT volume which matches a collagen fiber given both the OCT image features and the detected orientations acquired using our 3-D mapping method. This technique was previously demonstrated on en-face OCT images of human atrial tissue [12] and has been adapted to work for 3-D volumes for the purposes of this manuscript. We chose a particle step size of 15 pixels and generated 80 particles per step. The particle start points were evenly distributed in the en-face plane and placed just below the tissue surface in the axial dimension. The interested reader can refer to Gan, et al.’s manuscript for more details on the particle filter and its implementation.

3. Sample preparation and imaging protocol

3.1 Tissue collection

Five human uterine tissue samples, three non-pregnant (NP) and two pregnant (PG), were collected from consenting patients following hysterectomy according to protocol approved by the Columbia University Institutional Review board. Samples were divided by cutting blocks of tissue at the anterior, posterior, and fundus. The blocks were sliced into 15-25 mm x 15-25 mm square slabs with 2-4 mm thickness, resulting in anywhere from 4 to 8 slices per tissue block. Tissue was flash frozen using dry ice and stored in a −80 C freezer. Because hysterectomy is performed in response to disease of myometrial tissue, our specimens represent a rich diversity of pathologies and parity to provide a heterogeneous dataset for testing our algorithm. Patient demographic information can be found in Table 1.

Table 1. Patient specific obstetric history^a

View Table

3.2 OCT imaging

Prior to imaging, samples were thawed and submerged in PBS to keep the tissue hydrated. All samples were imaged using a commercial TELESTO SD-OCT system (Thorlabs, GmbH, Germany) with 6.5 $\mu$m axial and 15 $\mu$m lateral resolution. Multiple volumetric scans were needed to image each tissue slice. Each image volume was 5.5 x 5.5 x 2.51 mm (1375 x 1375 x 512 voxels) and overlapped 0.5 mm with adjacent volumes for digital stitching. Tissue slices were imaged front and back to capture the full sample thickness. The number of sub-volumes acquired varied by tissue slice depended on size and shape and ranged from 6 to 20 sub-volumes for each side of a tissue slice. The wide range in the number sub-volumes is a result of the varying tissue block’s size. Sub-volumes were stitched using a custom OCT registration and blending method to create a single mosaic volume [32]. This resulted in 28-41 mosaic volumes per patient.

4. Results

4.1 3-D fiber simulation

The proposed 3-D fiber orientation method was first demonstrated on a simulated image volume. Collagen fibers were simulated by generating high intensity tracts within a noisy image volume. All simulated fibers were oriented at $\theta = 20^{\circ }$ and $\phi = 70^{\circ }$ and uniformly placed in the image volume with 10 pixel lateral spacing. Figure 6(a) shows a 3-D view of the simulated fibers. The fibers within the image volume are visualized by looking at orthogonal image planes which are shown for the en-face plane (Fig. 6(b)) and the transverse plane (Fig. 6(c)) and were taken at the locations represented by the red and green planes in Fig. 6(a), respectively. Figures 6(d) and 6(e) show the fiber orientation maps for $\theta$ and $\phi$, respectively. The uniform color of the maps indicate the uniform distribution of fiber orientations. Some error in measurement was observed largely due to effects from the discrete nature of the thin, simulated fibers. This resulted in an average error of $2.58^{\circ }$ in $\theta$ and $0.04^{\circ }$ in $\phi$.

 

Fig. 6. 3-D fiber orientation mapping example for simulation image volume. (a) Simulated fibers generated at $\theta = 20^{\circ }$ and $\phi = 70^{\circ }$ and evenly distributed throughout an image volume. (b) En-Face (XY) view of image volume (red plane in (a)). (c) Side view (XZ) of image volume (green plane in (a)). (d) Map of the depth-averaged $\theta$ orientation measured using the proposed fiber orientation algorithm. Total average orientation measured was $22.58^{\circ }$ ($2.58^{\circ }$ error). (e) Map of the depth-averaged $\phi$ orientation measured using the proposed fiber orientation algorithm. Total average orientation measured was $69.96^{\circ }$ degrees ($0.04^{\circ }$ error).

Download Full Size | PDF

In biological tissue, collagen fibers rarely group in perfectly ordered bundles as they are in Fig. 6. Using a similar simulation, the 3-D mapping algorithm’s accuracy was assessed for fibers with different degrees of variation in their orientation. Fibers were simulated at random orientations with mean angles of $\mu _\theta = 20^\circ$ and $\mu _\phi = 70^\circ$, and some standard deviation $\sigma _\theta$ and $\sigma _\phi$. $\sigma _\theta$ was varied from $0^\circ$ – $90^\circ$ while $\sigma _\phi$ ($\sigma _\phi$) was varied from $0^\circ$ – $45^\circ$. For each simulated pair $(\sigma _\theta ,\sigma _\phi )$, the absolute value error in orientation was calculated for $\theta$ and $\phi$ ($e_\theta$, $e_\phi$) and averaged over all voxels of the simulated image volume using circular statistics. Figure 7(a), 7(b), 7(c) shows the simulated fibers for three different $(\sigma _\theta ,\sigma _\phi$), pairs. As $\sigma _\theta$ and $\sigma _\phi$ increase, the orientation of adjacent fibers becomes more varied and the full fiber group becomes more chaotic. Figures 7(d) and 7(e) show the mean absolute value error for $\theta$ and $\phi$, respectively. For larger values of $\sigma _\theta$ and $\sigma _\phi$ the error is higher due to increasing number of overlapping fibers and variance in orientation within algorithm’s measurement windows. The maximum value of $e_\phi$ was $31.05^\circ$ for $\sigma _\theta = 70^\circ$ and $\sigma _\phi = 40^\circ$ while the maximum $e_\theta$ was $20.27^\circ$ for $\sigma _\theta = 90^\circ$ and $\sigma _\phi = 40^\circ$.

 

Fig. 7. Fiber proximity simulation mapping test. A group of uniform randomly oriented fibers were simulated with mean angles $\mu _\theta = 20^\circ$ and $\mu _\phi = 20^\circ$ and standard of deviation $\sigma _\theta$ and $\sigma _\phi$ (a) Simulated fiber bundle with all fibers at the same orientation of $\theta = 20^\circ$ and $\phi = 70^\circ$. (b) Simulated fiber bundle with moderate deviation in orientation from the simulation in (a) ($\sigma _\phi = 20^\circ$, $\sigma _\theta = 30^\circ$). (c) Simulated fiber bundle with heavy deviation in orientation from the simulation in (a) ($\sigma _\phi = 45^\circ$, $\sigma _\theta = 90^\circ$). (d) Absolute value error in $\theta$, $e_\theta$. (e) Absolute value error in $\phi$, $e_\phi$. $e_\theta$ and $e_\phi$ both increase with $\sigma _\theta$ and $\sigma _\phi$. $e_\theta$ is nearly independent of $\sigma _\phi$, while $e_\phi$ is affected most significantly by $\sigma _\theta$.

Download Full Size | PDF

The en-face orientation $\theta$ is calculated by analyzing local image gradients within a window of pre-defined size. It is known for window-based methods that their accuracy is dependent on the window size. This relationship was evaluated on simulated fibers by varying the window size for different degrees of fiber variance, similar to the test in Fig. 7. Again, simulated fibers were drawn at random orientations with a mean orientation of $\mu _\theta = 20^\circ$, $\mu _\phi = 70^\circ$ with three different standard of deviations (see Fig. 7 Legend). Window size was varied from 11 x 11 to 71 x 71 pixels in steps of 10 pixels and the average orientation ($\mu _\theta$, $\mu _\phi$) and error ($e_\theta$, $e_\phi$) were calculated for each window size. The results are shown in Fig. 8, where the average measured orientations $\mu _\theta$ and $\mu _\phi$ are plotted in Fig. 8(a) and Fig. 8(b) and the average errors $e_\theta$ and $e_\phi$ are plotted in Fig. 8(c) and Fig. 8(d). $\mu _\theta$ and $\mu _\phi$ as a function of window size was found to be dependent on $\sigma _\theta$ and $\sigma _\phi$. For the three simulation cases, $\mu _\theta$ and $\mu _\phi$ remained close to the true mean orientations ($\mu _{true,\theta } = 20^\circ$, $\mu _{true,\phi } = 70^\circ$), but the variance as a function of window size was higher for larger values of $\sigma _\theta$ and $\sigma _\phi$. Similarly, $e_\theta$ and $e_\phi$ remained relatively constant for small $\sigma _\theta$ and $\sigma _\phi$, but increased as a function of window size for larger $\sigma _\theta$ and $\sigma _\phi$. The results show that when all the fibers of a group are at a similar orientation, large or small window sizes will capture the same orientation on average with high accuracy. When the deviation in orientation between adjacent fibers becomes more drastic, larger window sizes contain more fibers at different orientations which results in higher error measurements, but similar average orientation.

 

Fig. 8. Length scale simulation test. Average orientation ($\mu _\theta$, $\mu _\phi$) and error ($e_\theta$, $e_\phi$) were measured for different window sizes and variance in local fiber orientation ($\sigma _\theta$, $\sigma _\phi$). The legend shows the three simulation cases that were tested. For large $\sigma _\theta$ and $\sigma _\phi$, $e_\theta$ and $e_\phi$ increase with window size (yellow), but remain constant with window size for small $\sigma _\theta$ and $\sigma _\phi$ (blue). $\mu _\theta$ and $\mu _\phi$ remain within $+/-5^\circ$ for all three simulation cases and window sizes greater than 21 x 21 pixels. (a) Average angle $\mu _\theta$. (b) Average error $e_\theta$. (c) Average angle $\mu _\phi$. (d) Average error $e_\phi$.

Download Full Size | PDF

4.2 Image enhancement

To visualize and identify fibrous features in real OCT images rather than simulated image volumes we developed a custom fiber enhance pipeline. Each OCT image volume was pre-processed using this image enhancement pipeline prior to fiber orientation analysis. The visual effect of the pre-processing is shown in Fig. 9 for an OCT b-scan from an NP (Fig. 9(a)–9(b)) and a PG uterine specimen (Fig. 9(c)–9(d)). The NP example is from the anterior location while the PG example is from the fundus. Both examples were taken from a myometrial tissue slice. The raw b-scan in Fig. 9(a) was intensity adjusted for visualization while Fig. 9(b) shows the same image when fully processed. While fiber architecture is visible in the raw image, it is difficult to resolve individual fibers which are obscured by speckle, low contrast, and reflection artifacts (bright vertical "streaks"). In the processed image these problems have been mitigated, improving the visibility of collagen fibers. The image enhancement has a similar effect on the PG example where multiple collagen fiber groups can be visualized which are nearly impossible to discern in the raw image.

 

Fig. 9. Example of image enhancement pre-processing pipeline. (a) B-scan from a NP specimen (anterior, myometrium). (b) Enhanced version of (a). (c) B-scan from a PG specimen (fundus, myometrium). (d) Enhanced version of (c). Scalebar = 0.5 cm.

Download Full Size | PDF

4.3 3-D fiber modeling

3-D fiber analysis was initially applied to a 2.85 mm x 2.85 mm x 2.51 mm sub-volume from an NP uterus specimen. The sub-volume was taken from a region where fibers had relatively uniform alignment so that the mapping and tractography results could be visually compared with ease. 3-D views of the OCT sub-volume and fiber tractograph are shown in Figs. 10(a) and 10(b), respectively. Figure 10(c) and 10(f) show averaged views of the OCT sub-volume from two different directions. The resulting depth-averaged $\theta$ and $\phi$ maps are shown overlayed onto the XY view in Fig. 10(c) and Fig. 10(f), respectively. Corresponding XY and XZ views of the 3-D particle filter fiber model are shown in Fig. 10(e) and Fig. 10(h), respectively. Comparing the results of both the mapping and particle filter to the OCT, the general direction of the fiber tracts align with the fibers in the OCT image. Visualization 1 shows the fiber tracts overlayed onto the OCT sub-volume.

 

Fig. 10. 3-D fiber model validation using OCT uterus image sub-volume. (a) 3-D view of the OCT image sub-volume spanning 2.85 x 2.85 x 1.5 mm volume. (b) 3-D view of fiber tracts in (a) using particle filtering technique. (c) OCT sum-volume projection in the XY (en-face) plane. (d) Depth-averaged en-face ($\theta$) fiber orientation map. (e) En-Face (XY) view of the 3-D particle filter model. (f) OCT sum-volume projection in the XZ (transverse) plane. (g) Depth-averaged incline($\phi$) fiber orientation map. (h) Side (XZ) view of the 3-D particle filter model. Fiber tract colors were chosen randomly to enhance visual contrast. Scale bar = 0.5 mm.

Download Full Size | PDF

3-D fiber orientation maps were obtained for two OCT mosaic volumes. Figure 11 shows the 3-D maps produced for a tissue slice from an NP (a,c,e,g) and PG (b,d,f,h) specimen. The orientation maps are overlayed with transparency onto the OCT sum-volume projection to provide a visual connection to the imaged tissue features. The colors indicate the local fiber orientation for $\theta$ and $\phi$ according to the colormaps shown on the right side of the figure. The 2D maps (Figs. 11(a), 11(b), 11(e), 11(f)) are depth-averaged, meaning the displayed orientation is the circular mean of all measured orientations at a given en-face location and over the axial range of the OCT image volume. Additionally, the 3-D fiber orientation maps are shown in Figs. 11(c), 11(d), 11(g), and 11(h). Their shape matches that of the corresponding tissue slice. The two tissue specimens displayed in the figure are from approximately the same anatomical location, however, the fiber organization differs significantly between them. In both the NP and PG examples (Fig. 11(a)–Fig. 11(d)), the $\theta$ fiber map shows how the the collagen fibers tend to group locally, but also swirl and frequently change direction across the full span of the tissue. Conversely, the NP $\phi$ fiber map in Fig. 11(e) shows a uniform fiber alignment parallel to the tissue surface with a few areas of moderate incline. The PG $\phi$ map in Fig. 11(f) shows that this fiber network has a slightly higher incline on average than in the NP example. Additionally, there are more areas of angled or vertically aligned fibers in the PG sample.

 

Fig. 11. En-Face ($\theta$) and incline ($\phi$) collagen fiber orientation maps from a non-pregnant (NP) and pregnant (PG) uterine sample. There are two different visualizations of these maps shown: 2-D depth-averaged orientation maps which are overlayed onto the OCT sum-volume projection for visualization (a,b,e,f) and full 3-D block visualizations (c,d,g,h). The different colors indicate the collagen fiber’s en-face orientation (a,c,e,g) or its incline with respect to the tissue surface (b,d,f,h) according to the displayed colormaps where $\theta ,\phi \in [0^{\circ },180^{\circ }]$. Scalebar = 0.5 cm.

Download Full Size | PDF

Histogram analysis can be applied to the $\theta$ and $\phi$ fiber maps to obtain more information about the tissue’s collagen organization. To demonstrate this, we looked at histograms of the 3-D fiber orientation for an NP and PG uterine specimen with a highly organized collagen networks. Each histogram was generated using a bin width of $0.5^{\circ }$ and covering a range of $180^{\circ }$ ($\theta \in [0^{\circ },180^{\circ }]$, $\phi \in [-90^{\circ },90^{\circ }])$. The NP sample histograms of the $\theta$ and $\phi$ orientation for the full mosaic volume are shown in Fig. 12(a) and 12(b), respectively. The $\theta$ histogram has three dominant peaks ($40^{\circ }$, $125^{\circ }$, and $175^{\circ }$), suggesting this volume has three main fiber groups, each defined by a dominant orientation and dispersion. These same groups, however, do not appear in the $\phi$ histogram in Fig. 12(b). Instead, it shows a single distribution centered at $0^{\circ }$ with the majority of non-zero bins in the range of $+/-40^{\circ }$. To further investigate the fiber groups, we examined the same histograms for three different sections of the volume. Figure 12(c) shows the $\theta$ map overlayed onto the OCT sum-volume projection and the three regions that were analyzed. These areas were roughly chosen to cover three areas with distinctly different appearance in the OCT en-face images. The resulting $\theta$ and $\phi$ histograms for each section are shown in Fig. 12(d). Sections 1 and 2 each cover a single dominant fiber group as indicated by the uni-modal appearance of the $\theta$ histogram and Section 3 has a mixture of different fiber groups. While all three sections have peaks in the $\phi$ histogram at $0^{\circ }$, each group has a different spread of orientations. Comparing Sections 1 and 2, the mean angle of Section 1 is lower and the variance is larger than Section 2. This suggests that the two groups may converge in the volume with the fibers from Section 1 weaving between and underneath the fibers in Section 2. The significantly higher variance of Section 3 compared with the other two sections suggests that multiple fiber groups may converge at that location, resulting in more a more chaotic fiber organization.

 

Fig. 12. Histogram analysis reveals fiber groups in the non-pregnant uterus sample. (a,b) Histogram of $\theta$ and $\phi$ orientations for full mosaic OCT image volume of an NP (anterior) uterus tissue specimen. (c) En-Face sum-volume image of the uterus specimen. Color overlay shows the depth-averaged $\theta$ orientation map. (d) Histogram analysis was divided into three separate sub-volumes, as indicated by the color outlines in (c). Each section covers a dominant fiber group with the corresponding histograms revealing the dominant orientation and dispersion of each group. Scale bar = 0.5 cm.

Download Full Size | PDF

The same analysis was repeated for a posterior tissue slice from a PG uterine sample and is shown in Fig. 13. The full volume $\theta$ histogram in Fig. 13(a) reveals up to four fiber groups with a more even spread between groups than in the NP example from Fig. 12. This is visually reflected in the $\theta$ map in Fig. 13(c). Again the volume was separated into three sections by visual identification of the fiber groups and histogram analysis was applied to each section (Fig. 13(d)). The $\theta$ section-based histograms have more exaggerated peaks for the different fiber groups, but overall, the sectioning does not isolate the different groups as well as in the NP example. Another interpretation of this result would be that the different fiber groups of this sample are more mixed and the overall organization is more chaotic than in the NP example. The $\phi$ histogram for the full volume in Fig. 13(b) again has no discernible grouping like in the NP example. When comparing the $\phi$ histograms for Section 1 to Sections 2 and 3 however, we see that the fiber incline characteristics change from the top half of the tissue slice to the bottom half.

 

Fig. 13. Histogram analysis reveals fiber groups in the pregnant uterus sample. (a,b) Histogram of $\theta$ and $\phi$ orientations for full mosaic OCT image volume of an PG uterus tissue specimen. (c) En-Face sum-volume image of the uterus specimen. Color overlay shows the depth-averaged $\theta$ orientation map. (d) Histogram analysis was divided into three separate sub-volumes, as indicated by the color outlines in (c). Each section covers a dominant fiber group with the corresponding histograms revealing the dominant orientation and dispersion of each group. Scale bar = 0.5 cm.

Download Full Size | PDF

Using the $\theta$ and $\phi$ fiber maps and the particle filtering technique, fiber tractographs were generated for the NP uterus sample shown in Fig. 12. Two different views of this tractograph are shown in Fig. 14. Fiber tracts were generated using a 3-D particle filter (green/blue) and overlayed over the 3-D OCT mosaic volume (gray). The different fiber tract colors distinguish the two largest fiber groups identified using the histogram analysis of the $\theta$ fiber map. The green fibers are characterized by $\theta \in [80^{\circ },155^{\circ }]$ while the blue fibers are characterized by $\theta \in [155^{\circ },30^{\circ }]$. Fibers tracts outside this range were omitted for visualization purposes. The two different views show how the different fiber groups interact. In the center of the volume, small patches of the green fiber group can be seen weaving between the blue group in a characteristic basket-weaving pattern. Animations of the fiber tracts alone are provided in Visualization 2 and Visualization 3.

 

Fig. 14. 3-D fiber model using particle filter technique for an NP sample. Color represents two dominate fiber groups identified via histogram analysis (see Fig. 12). a) View 1. b) View 2. Scale bar = 4 mm.

Download Full Size | PDF

5. Discussion

In this manuscript we demonstrated a complete imaging and processing workflow for analyzing and modeling 3-D collagen fiber architecture from SD-OCT volumes of human uterine specimens. Image processing pipelines were developed to enhance OCT image quality, acquire 3-D fiber orientation maps, and produce tractographic models of the uterine fiber network. The proposed method was demonstrated for simulated image volumes and ex-vivo specimens of pregnant and non-pregnant uterine tissue from multiple anatomical locations. We found that the proposed methods provide an efficient and automated workflow for both quantitative analysis and visualization of 3-D collagen networks.

Some similar methods have been previously proposed, most notably the method by Gan, to quantify myofiber direction in 3-D for SD-OCT image volumes of human heart specimens [12]. Our method differs from Gan, et al.’s method in a fundamental way. Their algorithm was specifically proposed for analyzing myocardial tissue where it is known that myofibers are highly organized and align parallel to the tissue surface with minimal deviation [33]. This assumption means that analysis can be restricted to en-face images and it was not necessary for them to measure fiber incline orientation directly from b-scans. While many tissues share this property, such an assumption is inappropriate for human uterine tissue where we found interweaving and vertically aligned collagen fibers to be prevalent. This may prove to be especially relevant in pregnant myometrium.

Other studies modeling 3-D fiber architecture in OCT have relied on the polarization information provided by PS-OCT to create fiber tractography models [21]. The techniques used in PS-OCT are powerful as they use both OCT structural information and polarization information to trace fiber tracts. Our goal in this study, however, was to explore methods for analyzing and tracing 3-D fiber networks from only the structural OCT information obtained using a common OCT system. We believe that by accomplishing this goal, 3-D fiber analysis can become more widely available to interested researchers who would like to use OCT systems that are less expensive, easier to use, and require less digital storage then PS-OCT systems.

A number of tests were used to validate the proposed method on simulated data. The tests investigated the accuracy of the 3-D orientation mapping for environments where fiber orientation varied greatly between adjacent fibers and for different en-face window sizes. For fibers with significant changes in orientation we found that the fiber incline $\phi$ was typically less accurate than $\theta$. In the worst case, the error in $\phi$ was as large as $35^\circ$. The error in $\phi$ becomes large in this scenario because the patch size used to sample the transverse images is larger than the window size used for $\theta$. The algorithm assumes a single orientation in depth, but this assumption will not hold when collagen fibers become chaotic on the length scale of a single image patch, which results in potential significant errors. In future iterations, modification of the patch processing method may allow for more granular sampling in the axial direction to account for multiple layers and overlapping fibers. These changes are essential to future statistical analysis.

The focus of this study was to develop and test a framework for 3-D fiber orientation mapping using OCT image volumes. While in-depth analysis of the resulting maps and tractographs was outside the scope of this study, histogram analysis was employed to demonstrate our method’s potential to reveal new information about uterine collagen macro-structure. Identifying fiber grouping on the millimeter length scale and analyzing their individual organization had not been previously reported and we feel this method of analysis will be impactful for understanding normal structure function as well as pathological uterine physiology as it pertains to patient-specific collagen organization. Analysis of fiber groups can be further expanded through Gaussian mixtures models to automate the process and make accurate measurements of mean fiber orientation and dispersion.

One interesting feature we observed in the histogram analysis (Fig. 12) is that the $\phi$ histograms in Figs. 12(b) and 12(d) have very sharp peaks at the $0^{\circ }$ bin. It is difficult to say what the source of this peak could be, however, it is likely to be a result of imperfections in the calculation of $\phi$. Because $\phi$ is calculated as the arctangent of a ratio (Eq. (1)), an orientation of $0^{\circ } \in [-90,90]$ will generally result when the numerator of the ratio is much smaller than its denominator. Additionally, $\phi$ is a function of $\theta$ and $\phi _{x,y}$, so small errors in these measurements will propagate through the $\phi$ calculation and may concentrate near $0^{\circ }$. Using a very narrow bin width would smooth this erroneous peak, but would also make the rest of the histogram impossible to interpret. We found that a bin size of $0.5^{\circ }$ was the best compromise between reducing the size of the $0^{\circ }$ peak and producing a smooth and interpretable histogram.

Uterine tissue is complex because of its sophisticated collagen fiber architecture and heterogeneous composition of cells and ECM components.. A limitation of this study is that the fiber orientation mapping and tractography phases do not include tools to distinguish different tissue types. Uterine tissue is comprised primarily of collagen and smooth muscle cells (SMCs), though relative concentrations of each are not well characterized and may change significantly between specimens. SMCs can be further characterized as either uterine smooth muscle (USM) or vascular smooth muscle (VSM) which are phenotypically different and vary in their presentation based on the presence of vessels and arteries [5]. Histological sectioning suggests that in general SMCs are also aligned structures and well correlated with collagen alignment [28]. This heterogeneity had a minimal effect on the accuracy of our 3-D fiber processing method because both quantification tools, the Radon method and the gradient method, inherently detect brighter image features which correspond with the collagen fibers. Erroneous measurements will only occur in regions with vary high concentrations of SMCs relative to collagen. In these scenarios, texture analysis or computer vision tools could easily be incorporated into future versions of this method to create models that capture all the heterogeneous features of the myometrium including USMs, VSMs, vasculature, and even pathological tissue.

When considering the uterine tissue complexity and heterogeneity, it is important to address the diversity of the samples used in this study. All uterine samples obtained and imaged using OCT were obtained from donors with different parity following hysterectomies. This implies that these samples are very likely pathological and do not represent the tissue features of a healthy individual. In preliminary and qualitative analysis of the OCT uterine images, we found that b-scans revealed a number of basic phenotypes regarding fiber incline. Some samples appeared highly ordered with fibers primarily at horizontal or vertical inclines while others had more chaotic organization with fibers oriented at seemingly random orientations. Visualization 4, Visualization 5, Visualization 6, Visualization 7, and Visualization 8 are videos showing a few examples of these phenotypes. It is still unclear how representative these images are of typical uterine tissue and further investigation is needed to more accurately assess the typical structure of the uterine ECM.

Despite these limitations, the proposed method offers significant advantages over previous methodologies for modeling the 3-D architecture of human uterine tissue. Lutton, et al. used histological processing to model a single NP human uterus 7 cm x 3 cm x 0.35 cm tissue block and create the highest resolution model produced over a multi-centimeter length scale [28]. OCT provides micron resolution images similar to histology and requires significantly less resources to image the same volume of tissue. Lutton required thousands of histology slides and an intensive registration method to create their model. Using OCT, we are able to image and fully model a tissue slice in a few hours. In combination with processing methods like our proposed fiber modeling approach, OCT holds promise for even more sophisticated structure function analysis. We envision that the minimal resources needed to employ our method will open the door to multi-cohort studies that would have been prohibitively expensive to accomplish using histological processing or similar methods. The interested reader will be able to download the 3-D fiber mapping code from Columbia University's Academic Commons.

6. Conclusion

An image processing method was developed for extracting and visualizing 3-D fiber orientation information from OCT images. The resulting fiber model is a global view of whole tissue fiber networks in human uterine samples which provides a patient-specific platform for examining collagen micro-structure on a macro-scale and comparing the modeled architecture between pregnant and non-pregnant specimens. The proposed method has the potential to provide the first micron resolution model of a full human uterus sample and valuable insight on the physiological mechanisms related to structural remodeling during pregnancy.