1. Introduction

During pregnancy, the uterine cervix is a mechanical barrier, retaining the fetus until term, which is 37 weeks of gestation. In some cases, the accelerated remodeling of cervical tissue and subsequent reduction in the mechanical integrity of the cervix have been hypothesized to be associated with preterm birth [1, 2]. The cervix is a dense collagenous tissue. The mechanical integrity of the tissue is in part attributable to the preferred directionality of its collagen fibers [3]. Mechanical properties of human cervical tissue have been previously measured using standard uni-axial compress/tension tests [4, 5], indentation [6], permeability tests [7], and aspiration [8–10]. During pregnancy, the cervix becomes softer to accommodate delivery [9]. It is hypothesized that the remodeling of the cervical collagens is the main feature of this softening process. In animal models of pregnancy it has been shown that this softening is accompanied by a decrease in collagen crosslinking [11]. In human tissue it has been shown that pregnant tissue is more soluble in weak acid solutions compared to nonpregnant tissue [12]. Accompanying this collagen crosslink turnover, it is hypothesized that the overall directionality (i.e. ultrastructure) of the cervical collagen becomes more disorganized as pregnancy progresses.

To understand the role of the cervical collagen ultrastructure during pregnancy, and as a complement to our mechanical testing and the collagen crosslinking studies, we present an optical imaging method to provide the collagen ultrastructural information of the cervix. There have been many imaging modalities [13–16] used to study the collagen structure of the cervix. Cervical structure from 8 nonpregnant human cervixes was initially investigated by Aspden, who proposed a three-zone model based on the measurements from X-ray diffraction [13]. One of the main features of the three-zone model [13] is the middle zone of circumferential collagen fibers that trace around the inner canal (e.g. the cervical canal that the fetus passes through during delivery). Additionally, more recent magnetic resonance image (MRI) [14] and second harmonic generation (SHG) data also describe a region of circumferential collagen fibers in the nonpregnant cervix [16, 17]. Currently, the size of the region of these circumferential collagens and the amount of fiber dispersion remain unknown, and the collagen organization of the pregnant human cervix is completely unknown.

Optical coherence tomography (OCT) has been previously used to image the cervix for cancer detection and staging [18–20] by analyzing changes of the layers within the cervix, including the epithelium, basement membrane, and stroma. The analysis of cervical ultrastructure, however, is still lacking, largely due to a major challenge that the diameter of a normal cervix (~30 mm) is larger than the field of view (FOV) of standard OCT systems. To achieve the objective of characterizing axial slices of the human cervix ex vivo, we will collect multiple OCT volumes and combine them through image processing.

Characterizing large samples with OCT imaging can be achieved by stitching multiple overlapping 3D volumes. Many of the current methods for stitching OCT image sets have been developed for retinal imaging, where blood vessels with high contrast are a good reference for registration. In addition, retinal tissue has a clearly layered structure, which facilitates registering B-scans in 2D. It is thus feasible to use SIFT to extract features from OCT images for registration [21]. Zawadzki et al [22] registered multiple retinal OCT volumes through manual stitching, then developed a ray cast method to render the volume and used an annealing algorithm to optimize the alignment among multiple retinal imaging volumes [23]. Multiple OCT volumes are stitched in [24] based on a blood vessel matching algorithm. In [25], a B-spline based free form deformation method was used to register multiple volumes for a motion-free composite image of the retinal vessels. In comparison, the study in [26] has stitched bladder tissues, which lack high contrast structures. The bladder sample was stretched to mimic normal distention. However, in many biological samples such as cervical tissue, the surface topology varies greatly increase the difficulty of stitching. To reduce the effects caused by varying topologies, we are building on the work for registering OCT images with low contrast features by developing an auto error correction method in the offset estimation between volumes and by employing gain compensation and multiband blending post processing steps.

In this paper, we present an automated registration method to stitch multiple OCT volumes to enable the study of the ultrastructure and collagen fiber dispersion of human cervical samples. We investigate the cervical ultrastructure, fiber directionality, and dispersion from pregnant (PG) and non-pregnant (NP) patients.

2. Methods and materials

2.1 Tissue collection

All samples were collected from an Institutional Review Board (IRB) approved protocol at the Columbia University Medical Center, where three PG and ten NP patients were consented. Three to four mm thick whole axial slices of cervical tissue, cut perpendicular to the inner canal, were obtained after the patients underwent a hysterectomy. The second slice from the internal os of the cervix was used in this study, where other tissue slices will be analyzed in future studies. The average ages of PG patients and NP patients were 29 and 42.7, respectively. Hysterectomies were performed on NP patients for benign indications and on PG patients for suspected abnormal placentation (accreta/increta/percreta). These cesarean hysterectomies were performed prior to the onset of labor. Samples were analyzed after clearance by a Pathologist.

2.2 Image protocol

Three-dimensional volumetric image-sets were obtained from human cervical samples imaged ex vivo. OCT image-sets were acquired using a commercial spectral domain OCT system, Telesto (Thorlabs GmbH, Germany). The system is an InGaAs based system with a source centered at 1325 nm and a bandwidth of 150 nm. The axial and lateral resolutions are 7.5 μm and 15 μm in air, respectively. The maximum axial line rate is 92 kHz. In our experiments, each volume consists of 800 × 800 × 512 voxels, corresponding to a tissue volume of 4 mm × 4 mm × 2.51 mm. Samples were placed on a linear translation stage underneath the objective. For each sample, we obtained multiple volumes. The sample was moved along the x-axis or the y-axis on the translation stage. Since the surfaces of the sample were not flat, the axial position of the sample was also adjusted to make sure that the sample was well focused. Therefore, there were offsets in the x, y, and z directions. There was an overlap of 10% to 20% between two adjacent volumes. Each volume has a corresponding color-camera image, delineating the FOV on the tissue. The camera images and OCT images were taken simultaneously. The two images were calibrated using the factory settings of the commercial system. In the calibration, a scaling factor is measured between camera image and OCT image in the x-y plane.

2.3 Algorithm flow of registration

Stitching of OCT images is an ongoing engineering challenge. Conventional registration methods such as salient feature region (SFR) [27] and speed up robust feature (SURF) [28] are not sufficient because the contrast within OCT images of biological tissue is much weaker than the contrast within other image modalities such as computed tomography and MRI. Though imaged at the same area, the extracted features from multiple overlapped OCT images may vary due to noise and/or axial position of the sample. Attempts to register images using existing software like XuvTools [29] are likely to fail due to the low contrast in OCT images. In addition, the volumes should be aligned in three dimensions, which increases the difficulty of the global optimization problem.

To stitch multiple volumes, our algorithm consists of three steps, as shown in Fig. 1: 1) registration within the en face plane, 2) registration along the z-axis, and 3) a post processing procedure. In Step 1, the registration is based on camera images. In Step 2 and Step 3, we use OCT images for the registration and post processing. For the first step, a) we pair individual volumes based on their scale invariant features [30] in the camera images. We measure the offset between each pair of volumes. b) Then we formulated a linear regression model between the measured offsets and the centroid of each volume in a global coordinate system. The least square estimation of the regression model is the global optimization registration results within the en face plane. c) An auto-correction method is then used to eliminate the erroneous estimates in pairwise offsets. For the second step, based on the estimation of offsets within the x-axis and y-axis, a) we calculate the overlapped area between pairs of volumes and measure displacements along the z-axis by comparing the measured edge in each volume. b) Another linear regression model is used to obtain the global offsets along the z direction. Based on the global offsets in the x-, y-, and z- axes, all volumes are stitched and visualized in 3D. c) Error auto-correction is used to eliminate erroneous estimation along the z-axis. Lastly, for the third step, a post processing procedure including a) gain compensation and b) multiband blending is used to smooth the transition band between adjacent volumes.

 

Fig. 1 Flowchart of the registration algorithm. Scale invariant transform (SIFT) is used in pair matching within the en face plane. Calibration along z-axis is based on estimation in first block and edge detection results. Linear regression models are used and least square estimations are calculated in global registration within en face plane and along the z directions. An error auto-correction method was used to detect and eliminate the unreliable estimation in pairwise estimation. After all offset are calculated, we use gain compensation and multiband blending method to reduce the artifacts which are caused by the inconsistent intensity among multiple volumes in the overlapped area.

Download Full Size | PDF

2.4 Step 1: Registration within the en face plane

The camera image of each volume is considered as the reference for the registration within the en face plane. The SIFT algorithm described in [30] is applied to each image. Hundreds of keypoints are extracted in each image automatically. A Best-Bin-First (BBF) algorithm [31] is used to search for matched keypoints within other images. A normalized descriptor, DES, with a dimension of 128, is assigned to each keypoint. The local descriptor is based on the gradient in a small area and thereby was robust to illumination and brightness changes in the RGB camera image. DES_i and DES_j are the descriptors of keypoint sets in image i and image j, respectively. Given a descriptor of a keypoint des ∈ DES_i, the set of distance (Dis_des) is defined as,

 

Fig. 2 Schematics of registration (a) in three dimensions; (b) within the en face plane of the camera image; (c) overlapped area in 3D; (d) overlapped B-scans. Within en face plane, the offsets between two volumes are measured through keypoint matching in (b). Here, only keypoints that are matched are shown. There are hundreds of keypoints in the camera image. The number of extracted keypoints and ratio of matched features vary with camera images and overlapped areas. According to our experiment, only a small proportion (~5%) of keypoints can be matched. The centroids of multiple volumes are estimated globally by considering their offsets. To register the OCT volumes along the z-axis, the B-scans at the overlapped area are compared in (d). The displacement of the detected edge shows the offset along the z-axis among volumes.

Download Full Size | PDF

Next, we stitch all the volumes in a global space. Denoting the centroid of the volume i (i = 1, 2, …, N) as (x_i, y_i), we construct a vector c = [x₁ x₂ … x_N y₁ y₂ … y_N]^T, representing the centroids in all volumes. Within the en face plane, the sample is moved linearly along x-axis and/or y-axis. Therefore, the offsets between two adjacent volumes are maintained in a global space. The geometric relationship between offsets and centroids of OCT volumes is shown in Fig. 2(a). Suppose there are M pairs of volumes. The measured offset, d = [dx₁ dx₂ … dx_M dy₁ dy₂ … dy_M]^T, has a linear relationship with centroid vector c:

Though the images are paired based on SIFT features, there is a possibility that the estimated offset may be erroneous due to the mismatched SIFT features in poorly featured images. If most of the offsets between matched keypoints are spurious, a bad link will be constructed between two images. In prior work, Random sample consensus (RANSAC) algorithm [32] was used to eliminate the outlier (bad links) in the paired images. However, RANSAC requires prior knowledge of the outlier scale, which is difficult estimate in our poor featured images. Rather than relying on a probabilistic model of an outlier scale, we propose an error auto-correction method based on a deterministic metric to detect and eliminate the erroneous pairwise estimates. Suppose the centroid of each volume is a node in 2D space and the estimated offsets are the links connecting all nodes. Considering a link connecting node i and j, we use a metric e(i, j) to measure the mismatch between pairwise estimation e_p(i, j) and globalized estimation e_g(i, j). The metric is defined as:

Theoretically, if the estimation of e_p(i, j) is erroneous, the globalized estimation e_g(i, j) should be quite different due to error accumulation. Also, if estimation e_p(i, j) is accurate and the global optimization is computed based on reliable nodes, the difference between e_p(i, j) and e_g(i, j) will vanish.

The error correction method consists of two steps: a pruning process and a spanning process. The pruning process constructs a reliable set of nodes to ensure an accurate global estimation for each link by reducing unreliable nodes. The spanning process enlarges the number of nodes in the tree while keeping high accuracy, according to the Dijstra algorithm [33].

Algorithm 1 outlines our method to prune the original trees. In each iteration, we examine the error e(i, j) and delete the node that causes the largest error. The iteration ends when the maximum error is smaller than a threshold. This results in a reliable subset of nodes.

Algorithm 1. Pruning process. Starting from all nodes, the algorithm progressively eliminates the node that causes largest error. The nodes in the list of good nodes construct the most reliable tree structure with all pairwise links at high accuracy.

View Table | View all tables in this article

Algorithm 2 outlines our method to add nodes to the reliable subset. In each iteration we add a new node, which is from original data set. All newly created links are examined and the link that causes the large error will be deleted. The algorithm terminates upon addition of all nodes. During the spanning process, all reliable links are maintained.

Algorithm 2. Spanning process. Starting from the list of good nodes, the algorithm iteratively adds nodes and links from the list of bad nodes. When one node is added, every newly constructed link is examined. Links that cause large error are eliminated.

View Table | View all tables in this article

An illustration of error auto-correction is shown in Fig. 3. There are 53 nodes in the original tree at the beginning as shown in Fig. 3(a). The magnitude of the link shows the displacement. The colors of links refer to e(i, j) of the link. Blue represents an error smaller than 1 pixel while red represents an error larger than 1 pixel. At each step, the node causing maximum e(i, j) is deleted. The pruning process does not terminate until the maximum e(i, j) of the remaining links drops below a set threshold. Figure 3(a) to Fig. 3(c) shows the pruning process, where nodes in Fig. 3(c) form a set of nodes that includes all reliable nodes and links.

 

Fig. 3 Schematics of error auto-correction within the en face plane. Each node denotes an OCT volume within the en face plane. Each link represents a pairwise offset estimation based on SIFT. The links are color coded based on the error between pairwise offset estimation and global optimization, e (i,j). Blue color indicates that the error is within 1 pixel, and red indicates that the error is over 1 pixel. The pruning process is shown from (a) to (c). The initial structure is in (a). After a deletion of 8 nodes, only a small number of erroneous links exist in (b). The most reliable tree is established after deleting 17 nodes as shown in (c). The spanning process is shown from (c) to (e). Nodes are added at each iteration. Only reliable links are added to the tree structure at each iteration. Therefore, the structure of (d) is more reliable than that of (b), though the two figures share the same number of nodes. The constructed tree is shown in (e). With pruning and spanning, we delete the erroneous links and construct an offset tree with highly reliable links.

Download Full Size | PDF

Based on the reliable set of nodes, a spanning process is executed by adding new nodes and links, as shown in Fig. 3(d) to Fig. 3(e). At each step, one node is added to the set of nodes. The links between newly added nodes and existing nodes are examined based on the measurement of e(i, j). If the link causes large e(i, j), the link will be deleted from the tree structure.

2.5 Step 2: registration along the z-axis

Based on the estimation of $\widehat{c}$, we can arrange the volumes as overlapped tiles in a single space. However, the volumes are not perfectly stitched because there are variations along the z-axis due to adjustment of the sample along z direction during imaging. To register multiple volumes along the z direction, we use the surface edge within the B-scan as a reference. We calculate the overlapped region between volume pairs. A typical example is shown in Fig. 2(c). To calibrate the displacement along the z direction, we obtain B-scans from two volumes at the same position within the overlapped region at two orthogonal directions (Fig. 2(c)). The tissue surface edge within each B-scan is determined by searching for the maximum gradient after smoothing the image with a median filter. The displacements of two edges are measured along the x-axis or the y-axis. A typical scenario of overlaid edges is drawn in Fig. 2(d). The measured offset ∆z_xij between volume i and volume j is computed as the median value of all displacements between edges in volume i and j along the x-axis. Also, we can measure the ∆z_yij between volume i and volume j along the y-axis. If two estimated displacements agree with each other, e.g., within 5 pixels, we take the average of the two estimations as the final offset between volume i and j. Similarly to the registration within the en face plane, a linear relationship can be formulated between the z coordinates of the centroid in all volumes, z = [z₁ z₂ … z_N]^T, and the measured offset, d₂ = [dz₁ dz₂ … dz_M ]^T as following

Similarly, we process the error correction along the z-axis. The mismatch between pairwise estimation and global estimation is calculated. Pruning and spanning process are employed to minimize the aforementioned mismatch. Unlike the case within the en face plane, the mismatch is computed in one dimension, the z-axis.

2.6 Step 3: Post processing

With all offsets between volumes globally estimated, we next combine the volumes into a stitched volume. In the overlapped region among multiple volumes, the contrasts vary significantly, especially for samples that have distinctive topology. In that case, there will be a visible edge in the combined volume. To solve this issue, we use gain compensation and multi-band blending [30, 34] to normalize and smooth the combined volume. Gain compensation is to normalize the brightness of each volume and multi-band blending is to smooth the edge between multiple overlapped volumes. A schematic of the post processing procedure is depicted in Fig. 4.

 

Fig. 4 Schematics of post processing in two dimensions. For all input images, the overlapped region is calculated. Then a gain factor for each volume is computed based on the intensity of overlapped region. The intensity of each volume is tuned by the gain factor. Each B-scan is divided into multiple bands. The banded images are weighted with different matrix based on the idea of blending low frequency band with larger range and blending high frequency band with smaller range. The reconstructed image is the weighted summation of all bands.

Download Full Size | PDF

Gain compensation

We denote g_i, g_j as gain factors for the overlapped volume i and volume j. To estimate the gain factor, we define an error function of gain g based on the number of the overlapped voxels and their intensities as following [30]:

Multiband blending

After gain compensation, the volumes are normalized to achieve a uniform brightness. However, the edge of the combined volumes can still be distinguishable due to the difference in SNR levels within overlapped volumes, where the surfaces within the overlapped regions have different imaging depth. Here, we use a multiband blending technique to smooth the edges. The main idea of multiband blending is to divide the original image into multiple bands and to assign a different weight to each band. For each overlapped image, we assign a weight function, W(x,y) = w(x)w(y), where w(x) varies linearly from a value of 1 at the center to a value of 0 at the edge. The image then has a weight matrix as $W^n(i,j)$ defined as the following:

Both the original image and weight matrix are divided into multiple bands through convolution with a Gaussian filter by computing the following in Eq. (12)-(14):

2.7 Fiber orientation estimation

Using the above techniques to enlarge the field of view, we can investigate the collagen fiber orientation and dispersion within the whole cervical axial slices from NP and PG patients. The gradient based algorithm that was previously used to quantify fiber orientation in the myocardium was used to quantify collagen fiber orientation [35, 36]. After the stitched volume is obtained, we determine the three-dimensional surface of the whole volume. Then, we obtain a plane parallel to the surface of the volume. Ten images are averaged along the z direction to reduce noise and improve the measurement of collagen fiber orientation. The fiber orientation algorithm is then applied to the plane. Briefly, the fiber orientation algorithm is as follows. To sharpen the image, a second order Butterworth high pass filter is convolved with the averaged image. In addition, a median filter is used to reduce speckle noise. For each image pixel (i, j), the intensity gradients in the horizontal (G_x) and vertical (G_y) direction are calculated by convolving two 3 × 3 Sobel. The magnitude of gradient G(i, j) and angle Φ(i, j) are calculated as [37]:

3. Results

3.1 Performance of stitching algorithm

Using the same data set from Fig. 3, we demonstrate our auto-correction method by plotting the maximum error e(i, j) versus the number of nodes in the pruning process and the spanning process in Fig. 5. Of note, Fig. 5(a) is plotted using a decreasing number of nodes. We observe a drastic decrease in the maximum error with a decrease in the number of nodes. Within this example, when the number decreases from 53 to 52, a drastic change in error is observed because the node that is removed has a large number of the links that are inaccurate. If we include that node, most of the estimations are incorrect. If that node is deleted, the error is largely reduced. The error decreases to <1 pixel and a reliable set of nodes is obtained. In this example, the number of reliable nodes is 37. For the spanning process, the error stays at a low value because we only add reliable links for each new node within each iteration.

 

Fig. 5 Measured maximum error between pairwise offset estimation and global optimization during (a) pruning and (b) spanning process. The x-axis in (a) is in decreased numbers. The curve between 51 to 37 is zoomed in the inserted figure. In pruning, the measured maximum error drops with number of nodes in global estimation. Through link selection, the maximum error remains unaltered with additions of new nodes.

Download Full Size | PDF

We set up an evaluation test to validate the accuracy of our method. We obtained an OCT volume of 4 mm × 4 mm × 2.51 mm (800 × 800 × 512 voxels) of a plastic cap as ground truth. Then, we image the same space with four overlapped volumes of 2.5 mm × 2.5 mm × 2.51 mm (500 × 500 × 512 voxels). The volumes are stitched and compared with the ground truth. We set four pairs of landmarks at ground truth and then specify the same landmarks within the stitched volume. The 3D distances between pairs of landmarks are measured using Amira, as shown in Fig. 6. Difference in measured distance between stitched volume and ground truth is 11.64 ± 2.32 voxels, showing a good correspondence between stitched volume and ground truth.

 

Fig. 6 Comparison of (a) four stitched volumes and (b) the whole volumes. Four volumes with sizes of 2.5 mm × 2.5 mm × 2.51 mm were imaged and stitched using our algorithm in (a). Within the same space, another OCT volume with a size of 4 mm × 4 mm × 2.51 mm was acquired for comparison. Eight landmarks were specified in each volume, formatting four pairs of distances. The distances were measured and compared in voxel as a validation of the accuracy of stitching algorithm. The maximum mismatch in the two volumes does not exceed 15 voxels.

Download Full Size | PDF

The impacts of our registration algorithm are shown in Fig. 7, using results from a PG human cervical sample as an example. First, the impact of error correction is shown in panels a and b. The stitched en face image without error correction is shown in Fig. 7(a) and the stitched en face image after error correction is shown in Fig. 7(b). Without error correction, obvious errors are observed, where we can see discontinuity in the sample boundary in some areas. We have observed this in other data sets when the number of volumes to be stitched is large. Large errors exist and the contour of the sample is inaccurate. After the error correction step is applied, the offsets between volumes are well estimated and the contour of the cervix is very smooth, Fig. 7(b). A representative stitched B-scan is shown in Figs. 7(c)-7(e) to show the impact of the post processing step. Figure 7(c) shows the stitched B-scan with hard combination, where images are combined by averaging the overlapped regions. Figure 7(d) shows the stitched B-scan after gain compensation and Fig. 7(e) shows the results after both gain compensation and multiband blending. For the same area, marked as Fig. 7(c1), Fig. 7(d1), and Fig. 7(e1), the transition band from one B-scan to the other is very distinguishable in the results using hard combination. Generally, the pixel value on the right half of Fig. 7 (c1) is lower than the value on the left half. After gain compensation, as shown in Fig. 7(d1), the pixel values are generally decreased in the left half to match the pixel value in the right half. However, the boundary is still very clear. With both gain compensation and multiband blending, as shown in Fig. 7(e1), the transition between the left half and the right half is smooth and the boundary is nearly invisible.

 

Fig. 7 Registration results within the en face plane (a) with and (b) without the error correction. Registration of B-scan results from (c) hard combination of multiple B-scans, (d) gain compensation of multiple B-scans, and (e) gain compensation plus multiband blending of multiple B-scans. The volumes are not stitched properly without error correction in (a). After error correction, the volumes are well arranged within the en face place with a consistent contour of the cervix. In a zoomed area indicated by dashed yellow rectangle, if hard combined, the two sections approximating the edge have great difference in pixel value in (c1). Gain compensation uniformed the pixel value but a boundary is still visible in (d1). When proceeding with multiband blending, a smooth transition of pixel values is visible over the edge in (e1).

Download Full Size | PDF

The algorithm was applied to all 13 samples. Six representative three-dimensional reconstructions of the stitched volumes of cervical axial slices are shown in Fig. 8, along with the stitched camera image. Figures 8(a)-8(d) are of two cervical axial slices from two PG patients. In Figs. 8(a)-8(b), 53 volumes are stitched as a new volume of 30 mm × 27 mm × 7 mm. In Figs. 8(c)-8(d), 54 volumes are stitched as a new volume of 28 mm × 27 mm × 4 mm.

 

Fig. 8 Registered results of (a-d) PG and (e-l) NP samples. In particular, (a, c, e, f, i, k) are stitched camera image in two dimensions while (b, d, f, h, j, l) are stitched OCT volumes. For one PG sample, 53 volumes are stitched as a new volume of 30 mm × 27 mm × 7 mm. For another PG sample, 54 volumes are stitched as a new volume of 28 mm × 27 mm × 4 mm. We stitched 40 volumes for the first NP sample and formed a whole volume of 20 mm × 24 mm × 5 mm. We stitched 41 volumes for the second NP sample and formed a whole volume of 21 mm × 26 mm × 5 mm. We stitched 24 volumes for the third NP sample and formed a whole volume of 15 mm × 19 mm × 5 mm. We stitched 48 volumes for the fourth NP sample and formed a whole volume of 24 mm × 31 mm × 5 mm. Smooth surfaces are observed when multiple volumes are combined.

Download Full Size | PDF

Figures 8(e)-8(l) are of four cervical axial slices from four NP patients. In Figs. 8(e)-8(f), we stitched 40 volumes and formed a whole measuring 20 mm × 24 mm × 5 mm. In Figs. 8(g)-8(h), we stitched 41 volumes and formed a volume measuring 21 mm × 26 mm × 5 mm. In Figs. 8(i)-8(j), we stitched 24 volumes and formed a volume measuring 15 mm × 19 mm × 5 mm. In Figs. 8(k)-8(l), we stitched 48 volumes and formed a volume measuring 24 mm × 31 mm × 5 mm. The results demonstrate that our method can provide a full view of the cervical axial slice.

3.2 Analysis of fiber structure

The directionality maps are shown in Fig. 9 for both PG cervical samples (a-d) and NP cervical samples (e-l). For each sample, we determine the fiber orientation in each 1000 μm × 1000 μm sub-region. Within Fig. 9, representative OCT images are taken at 245 μm below the surface in all volumes. Results are shown at 10 depths below 245 μm with an increment of 49 μm in Media 1 and Media 2 for a typical PG sample and a typical NP sample. The color of the vector on fiber orientation maps were encoded with a confidence value defined Eq. (19). We observe that the fiber orientation in NP cervical samples are more regular with higher confidence. The fiber orientation in PG cervical samples are less regular and with lower confidence, especially at the region close to the inner canal. Moreover, we observed that the orientation does not vary within the depth between 245 μm and 735 μm regardless if the sample was PG or NP.

 

Fig. 9 Representative OCT images of (a-d) PG and (e-l) PG cervix. The FOV are 30 mm × 27 mm in (a)-(b), 28 mm × 27 mm in (c)-(d), 20 mm × 24 mm in (e-f), 21 mm × 26 mm in (g)-(h), 19 mm × 15 mm in (i)-(j), 24 mm × 31 mm in (k)-(l) . The image is a 2D plane that is parallel to the surface with a vertical depth of 245 μm. Fiber orientation results were processed for each OCT image. The estimations of orientation were made on each 1000 μm x 1000 μm. The results are color-coded based on confidence value. Two typical results are shown at 10 depths below 245 μm with an increment of 49 μm in Media 1 and Media 2 for a PG(a-b) sample and a NP(e-f) sample. The fiber orientation does not vary much within the range of depths.

Download Full Size | PDF

Figure 10 shows a dispersion analysis of collagen fibers within 1000μm × 1000μm × 490μm sub volumes (marked as boxes), corresponding to the approximate dimensions of individual elements used within Finite Element analysis for cervical ultra-structure [3]. The probability distribution of orientations calculated from Eq. (18) in twelve sub-regions over each axial slice is plotted. In particular, the regions located anterior, posterior, left, and right side of the inner canal are evaluated. Within each side, 3 locations are selected to represent inner, middle, and outer regions. As shown in Fig. 10, the inner region is defined as the first full window closest to inner canal. The outer region is defined as the last full window distant to the inner canal. The middle region is the window in the middle of inner region and outer region. Figures 10(a)-10(l) compared the distribution between a representative NP sample (red colored curve) and a representative PG sample (blue colored curve) at the twelve regions described above. The probability is plotted as a function of angle. In general, the distribution of NP sample has a higher peak value than the distribution of PG sample, which implies that the orientation of collagen fiber bundles show less dispersion and more regularity for NP sample. Generally, the sub-regions close to inner canal (k, e, f, and l) show larger dispersion than the distant sub-regions (a, b, g, and h).

 

Fig. 10 Distribution over 12 sub-regions (a) – (l) of non-pregnant (NP) (red) and pregnant (PG) (blue) samples. For each sub-region with a size of 1000 μm x 1000 μm area, the distribution is averaged over 10 depth with an increment of 49 μm. From (a) to (l), the peak value of distribution in NP is higher than the value in PG, which indicates less dispersion and more regularity.

Download Full Size | PDF

To further analyze the dispersion of fiber organization, the confidence value, which gives an indication of fiber dispersion, is evaluated at all 12 locations from all 13 cervical samples. A higher confidence value indicates lower dispersion. The mean and standard deviation are summarized in Table 1. We observe that the confidence from NP samples is higher than the confidence measured within the PG samples, indicating an increase in fiber dispersion within the PG cervical samples. In addition, a general trend is observed within both PG and NP samples, where the region adjacent to the inner canal has lower confidence (higher dispersion) than the middle and outer regions. We hypothesize that that the inner regions show larger dispersion because of the function of inner canal.

Table 1. Statistic of the confidence value at inner region (e, f, d, l), middle region (c, d, i, j), and outer region (a, b, g, h). The confidence value is measured as mean ± standard deviation.

View Table | View all tables in this article

4. Discussion

In this paper, we present a method to enable the study of the ultrastructure of the cervix with high resolution using OCT. We developed a stitching algorithm to increase the field of view to allow for imaging of whole axial human cervical slices. We then applied a fiber orientation algorithm to analyze the collagen fiber dispersion patterns within different areas of the cervix and between PG and NP samples. Together, we have developed a tool-set to analyze the ultrastructure of cervical fiber. We are able to visualize the circumferential orientation of collagen fibers within axial cervical slices. Specifically, we show initial evidence that the collagen network remodels in pregnancy, where pregnant cervical collagen fibers maintain an overall circumferential direction but are more dispersed when compared to the nonpregnant state.

Compared with the existing registration method used in retinal images [21–25], our work does not rely on any contrast structures. Our work can be readily applied to analyze other samples with low-contrast features. To register tissues lacking high contrast structures, efforts are made to stitch volumes in stretched bladder [26]. We have built upon this work by stitching large numbers of volumes within samples that were not stretched. Our correction method and multiband blending step help to ensure an accurate and smooth combination for large samples with varying surface topology. We showed that error correction is necessary especially when the number of volumes are large. Otherwise, there will be discontinuities in the images of the whole ultrastructure and it would be difficult to analyze the fiber orientation.

Our work has three limitations. First, the intensity decreases at higher axial depths. It is limited by the point spread function of OCT system along the axial direction. In a stitched en face image, the part imaged at a higher depth is less informative than the part imaged at a lower depth regarding to the morphological ultrastructure. Second, the samples were moved with linear translation stages. In the future, we will develop a generalized algorithm to take into account rotation, where a vector of translational displacement and a vector of rotation will be used to extract features to improve image set positions in the global space. Third, due to the limited availability of cervical samples from pregnant patients, we only collected 3 PG samples and there is an age difference between PG and NP samples in this study. For this reason, the statistical analysis had a lower power. Currently, we observe a trend that PG samples show larger dispersion compared to NP samples. This encourages further study to analyze if the dispersion is due to the pregnancy status of the patient or the age of the patient. With more cervical tissues samples, we will be able conduct a comprehensive analysis of the remodeling of cervical collagen during pregnancy, controlling for factors such as age, gestational week, and number of prior pregnancies.

Regardless of these limitations, the study of cervical collagen ultrastructure has been scarce, and OCT, along with our stitching algorithm with collagen fiber dispersion analysis provides a great tool ex vivo studies. In the future, we will use the directionality map and dispersion analysis of the whole axial slice as inputs into mechanical models of cervical tissue [38]. Importantly, the 1000μm sub-volumes used within our dispersion analysis correspond to the finite element sizes used within mechanical models of the cervix. Large scale fiber orientation maps, as shown in Fig. 9, specific to a tissue sample can also inform mechanical testing and provide a testbed for structure-function studies of whole axial cervical slices. It currently takes three hours to stitch forty volumes. The registration along z-axis and multiband blending account for most of the computational cost. In the future, we will implement real time stitching through distributed computing. With a reduction in the computational time, this framework can be used to analyze collagen fiber orientation and dispersion within axial slices of the whole cervix, where a typical sample will have 5-6 axial slices. This will allow us to investigate regional changes in the radial and axial directions, and with a larger sample size investigate factors that result in inter-patient variability. This will provide a wealth of information about the collagen ultrastructure of the cervix.

5. Conclusion

We observed increased collagen fiber dispersions within cervical samples from pregnant patients compared with cervical samples from non-pregnant patients. In addition, we observed regional differences in fiber dispersion in both PG and NP samples, where fiber dispersion was increased close to the inner canal. This was possible through an optimized image processing tool to investigate the cervical ultrastructure of collagen fiber using OCT. This image processing tool uses image stitching to enlarge FOV of OCT acquired 3D image sets. In addition, we mapped the directionality of fiber bundles of cervical axial slices and performed a dispersion analysis between NP and PG samples accordingly. Our work will enable future studies to quantify spatial variations in collagen fiber dispersion during the progression of pregnancy and analysis of fiber dispersion in other organs/tissues such as the myocardium. Furthermore, the view of cervical ultrastructure provided by enlarged FOV can potentially facilitate mechanical modelling of cervical samples with the aim of understanding preterm birth.