1.

Introduction

The extracellular collagen matrix (ECM) has been found to promote the progression of many types of cancer. However, the underlying mechanism behind this relationship is not fully understood, and is the subject of intense biomedical research. Much of this research has benefited from the capabilities of laser-scanning microscopy techniques, in particular second-harmonic generation (SHG) imaging¹ to capture high-resolution, high-contrast images of individual collagen fibers in tissue and in in vitro tissue models.^2–6 For example, Conklin et al.⁷ showed that patterns in SHG images of collagen can predict breast cancer patient outcome. Raub et al.⁸ showed that SHG image characteristics can be used to predict bulk mechanical properties of collagen hydrogels, a common in vitro tissue model for studying cancer cell motility. Nadiarnykh et al.⁹ and Watson et al.¹⁰ found that SHG image characteristics in ovarian tissue provide quantitative discrimination between tumor and benign tissues. Although SHG has been used successfully in these and many other studies, quantification of collagen fiber shape changes remains a difficult challenge, in part due to large heterogeneities in the patterns observed in SHG images of tissue. For example, in the breast tissue images shown in Fig. 1, collagen fibers can be described as wavy or straight [Figs. 1(a) and 1(b)], high or low density [Figs. 1(c) and 1(d)], and with thick bundles or thin strands [Figs. 1(e) and 1(f)]. These descriptions are consistent with previously published observations of common collagen structures in tissues.^7,11,12 In addition, depending on imaging parameters such as depth within the tissue, images can have low signal-to-noise ratio (SNR) and potentially low-dynamic range [Figs. 1(g)–1(j)]. Quantitative analysis techniques for SHG images of collagen need to provide robust and informative features within this heterogeneous collection of patterns and image qualities. Also, in order to elucidate the interactions between cells and individual collagen fibers, effective quantitative analysis techniques should be able to extract information about individual fibers such as fiber number, length, angle, and curvature. The work reported here is motivated by these two requirements: the need for robust performance and the need for fiber-level information in SHG image analysis of collagen.

Fig. 1

Representative collagen patterns observed in human breast cancer tissue sections demonstrating the heterogeneous nature of collagen structure. Wavy (a) and straight (b). High (c) and low (d) densities. Thick bundles (e) and thin strands (f). Discontinuous (g) and continuous (h). Crossing (i) and parallel (j). Scale bar is 10 μm.

While these two requirements may be met by manual analysis,¹³ inter- and intra-observer variance can be significant and time requirements can be prohibitive. Computer-assisted image feature extraction is poised to help meet these requirements in an automated fashion. Transform or filter-based methods have been used for SHG collagen analysis such as the Fourier transform method published by Falzon et al.,¹⁴ the combined Fourier and Hough transforms method by Bayan et al.,¹⁵ the curvelet transform (CT) method by Pehlke et al.,¹⁶ the Fourier and fractal-based methods reported by Frisch et al.,¹⁷ the directional gradient technique suggested by Altendorf et al.,¹⁸ and the gray-level co-occurrence method published by Hu et al.¹⁹ These techniques can be highly robust, often able to detect important features in a diversity of image settings. However, since these techniques do not extract individual fibers, they lack the ability to identify fiber-level information. For example, transform-based methods can provide general information about fiber size and direction at each point in an image but cannot determine the actual fiber number nor the length or curvature of the fibers. Since pixels are not grouped into individual fibers, these methods may not sense the difference between many randomly oriented, straight fibers and long curvy fibers, features that may help classify patients into high- and low-risk groups for ovarian⁹ and other cancers. In addition, angle distributions generated by these algorithms would generally produce bias toward longer and potentially thicker or brighter fibers, since more pixels are present in longer and thicker fibers and distributions are based on pixels and not fibers. These errors must be avoided in order to appropriately test the biological hypothesis that fiber angle distribution can help predict the metastatic potential of cancer cells.^4,5,7,20

Fortunately, fiber tracking and extraction methods, such as those published by Wu et al.^21,22 and Stein et al.,²³ have been developed to extract fiber-level information from images of in vitro collagen matrices. These methods can enable the automated measurement of important fiber-level parameters such as fiber length, number, and curvature, and have been used to estimate collagen gel mechanical properties based on confocal images of stained gels. However, they have not been applied to SHG images of collagen in situ. While these approaches are powerful, perhaps they have not been applied in situ because they often fail to properly segment fibers in the dense or low SNR situations commonly encountered in SHG images of tissue. Examples of two SHG images are shown in Figs. 2(a) and 2(d) with corresponding manual fiber extractions shown in Figs. 2(b) and 2(e). The fiber extraction (FIRE) algorithm, developed and made available by Stein et al.,²³ produces the overly complex fiber network shown in Fig. 2(c) and an erroneous star pattern in Fig. 2(f), and in both cases fails to identify many of the fibers extracted by the human observer.

Fig. 2

Fibers extracted by the FIRE algorithm alone without preprocessing. (a) and (d) are the original images, (b) and (e) show manual segmentations of the fibers, and (c) and (f) show the automatic fiber segmentations that are extracted by the FIRE algorithm and show many falsely segmented fibers. Scale bar is 25 μm.

Instead of using transform-based methods or fiber-extraction methods alone, a more strategic approach would be to combine complimentary methods and to use transform-/filter-based methods as a preprocessing step to fiber extraction. This combined approach has potential for robust performance in a wide range of challenging imaging situations commonly seen in cancer imaging, while simultaneously allowing for fiber-level information to be extracted from images. In this study, we present an approach for integrating transform-/filter-based preprocessing techniques with fiber extraction and evaluating the performance of our combined approach with the ultimate goal of improving the fiber extraction accuracy of an algorithm such as FIRE. Our hypothesis is that the application of appropriate pre- and post-image processing algorithms can significantly improve fiber extraction in SHG images. This will enable, for example, more accurate fiber angle distributions, thus allowing for increased sensitivity to detect collagen alignment changes related to cancer progression.

We evaluate four candidate preprocessing techniques including the simple Gaussian filter (GF), the SPIRAL-TV filter,²⁴ the tubeness filter,²⁵ and a CT-based denoising filter.^26,27 Other than the GF, these filters were chosen based on their published ability to highlight edge information in images, while simultaneously suppressing spatially uniform structures and noise. We have chosen to use the FIRE algorithm based on the evidence of its ability to extract fibers from in vitro collagen gel networks and its availability;²⁸ however, other fiber extraction tools may be substituted for the FIRE algorithm. We focus our analysis on two-dimensional (2-D) SHG images, because the effective nonlinear susceptibility declines sharply when fibers are tipped out of the imaging plane;^29–31 however, our methods may naturally extend to three-dimensional (3-D) without significant alteration.

2.

Materials and Methods

Our experimental approach for the evaluation of quantitative collagen fiber extraction is illustrated in Fig. 3. Twenty-five images were annotated by three human observers, who traced each fiber within each image to create a surrogate for ground truth. The same images are then filtered with one of four image filters, which are described in the following sections. Fiber tracking and extraction is then performed using the FIRE algorithm. We evaluate the automated fiber extraction accuracy by comparing the length, angle, and position of each fiber extracted by the FIRE algorithm with each fiber that was manually extracted to determine if fibers can be considered detected, missed, or falsely detected. $F$-measure scores are created based on these rates, and a single parameter in each image filter was adjusted to optimize the $F$-measure score for each algorithm. The details of each step in this process are given in the following sections. After algorithm evaluation, we applied the best performing technique to quantify collagen fiber shape changes during tumor progression in an in vivo mouse model for breast cancer.

Fig. 3

Diagram of the approach for quantitative collagen analysis showing the iterative process for optimizing the performance of a single image-processing filter for fiber tracking. The raw image is processed by the image filter using an initial normalization parameter, the result of which is sent to the FIRE fiber-tracking algorithm. Automated fiber extractions are compared against manually performed fiber extractions. Several normalization parameters are evaluated, and one optimal parameter is selected for each filter based on the fiber evaluation result.

3.

Results

Comparison of the four image-processing techniques to each other, as shown in row 1 of Fig. 4, reveals that edge-preserving filters, such as SPTV, although effective for denoising without the loss of edge information, do not lend themselves well to improving the fiber tracking results. On the other hand, the TF and CT create ridges along fiber centers (Fig. 4, row 1), helping to ease the difficulty of threshold selection and helping the fiber-tracking algorithm to follow the centers of thick or noisy fibers. Examination of fiber-tracking results in Fig. 4, row 2 shows many completely erroneous fiber tracks for the unprocessed, GF, and SPTV filtered cases (red arrows), whereas the TF and CT filtered results show several properly segmented fibers (green arrows). Each of the images in Fig. 4 is a representative 128 by 128 pixel region cropped out of larger images.

Fig. 4

Output of the image-processing techniques (row 1) and output of the fiber-tracking algorithm (row 2) for a single test case. The first column is without a filter, column 2: GF, column 3: SPTV filter, column 4: TF, and column 5: CT. Scale bar is 25 μm.

To further test the effect of preprocessing on collagen fiber segmentation, the manual segmentation results were compared with the automated fiber extraction results for two representative test cases, shown in Fig. 5. Each row in the figure is a different test case, while each column represents a different method of fiber segmentation. Column 1 shows the original images with no overlaid segmentations. Columns 2 through 6 show the original image with overlays of the manual, GF, SPTV filter, TF, and CT filter segmentations, respectively, where FIRE was performed following each of the filter preprocessing steps. Although we had three observers manually segment each of the test cases, the manual segmentations shown in column 2 represent the segmentations of a single observer. Each tile in Fig. 5 is a 128 by 128 pixel crop of a larger image. The quality of extraction of individual fibers under the various processing conditions was compared [arrows in Figs. 5(a) and 5(b)], where it was found that improperly segmented fibers as a result of GF or SPTV filtering were accurately segmented using the TF or CT filter.

Fig. 5

Two test cases (a and b), showing different processing methods in each column. The original image (column 1) is shown overlaid with a manual segmentation (column 2), GF (column 3), SPTV (column 4), TF (column 5), and CT filter (column 6), where each filter is followed by FIRE fiber extraction. Scale bar is 25 μm.

In order to assess the accuracy of collagen fiber segmentation, the results of each of the fiber extraction preprocessing algorithms were compared against each of the three segmentations performed by the independent observers using a collection of custom scripts written in MATLAB. If a fiber segmented by the automated process had a similar angle, close proximity, and similar length to a manually segmented fiber, then an association was made between the automated and manual fibers, indicating a true positive. After all fibers were evaluated, all remaining unassociated manual fibers were counted as false negatives (misses), and all remaining unassociated automatic fibers were counted as false positives (false hits). Precision, recall, and their harmonic sum ($F$-measure) were computed and compiled for all test cases and all observers. Overall average $F$-measure scores for each of the preprocessing algorithms are shown in Fig. 6. The average $F$-measure score for the CT filter was the highest followed by the TF, SPTV, and GF. The error bars indicate the standard deviation between the $F$-measure scores from each of the three observers and show that the scores between observers were very similar, meaning that the CT filter result was the closest match to all three observers. Comparing the computation time for each algorithm on identical images resulted in respective times of 0.22, 225.05, 0.83, and 4.65 s for the Gaussian, SPTV, Tubeness, and CT processing. Focusing on the CT denoising filter, we can further validate its performance combined with FIRE using a collection of computationally generated images of collagen fibers that were designed to mimic the length and curvature characteristics found in collagen gels of approximately $1.0\mathrm{mg}/\mathrm{mL}$.²³ We processed these images with the CT-FIRE algorithm to extract length and angle information about the individual fibers. The results of this test are shown in Fig. 7 and show that the CT-FIRE algorithm produces accurate length and angle distribution measures in known synthetic test cases.

Fig. 6

$F$-measure, recall, and precision results comparing the automated segmentation techniques with the manual segmentations of three independent raters, for 25 test cases, representing a total of 9290 fiber evaluations. The error bars indicate the standard deviation between average $F$-measure, recall, and precision scores of each of the raters. Recall is the fraction of relevant fibers that were found. Precision is the fraction of fibers found that was relevant. $F$-measure is the harmonic sum of recall and precision.

Fig. 7

Distribution of angles (top row) and lengths (bottom row) of all fibers in all simulated test cases. Ground-truth data is on the left, and the results of the automated CT-FIRE algorithm are shown on the right.

After identifying the CT filter as the top-performing preprocessing algorithm, we applied the combined CT-FIRE algorithm to the measurement of collagen fiber shape changes in an in vivo mouse model for breast cancer. The results of this study are shown in Fig. 8. Representative images show clear differences in waviness of fibers between the early [Fig. 8(a)] and late [Fig. 8(b)] time points. The colored lines overlaid on the images are the automated fiber segmentations produced by CT-FIRE. These overlays qualitatively illustrate the high-fiber segmentation quality that can be expected from the CT-FIRE algorithm. Fiber waviness ($W$) was quantified for each extracted fiber by dividing the distance along the path of the fiber ($d_0$) by the distance between the end points of the fiber ($d_n$).

Fig. 8

Demonstration of automated fiber segmentation feature extraction in an in vivo mouse model for breast cancer. A mammary window was placed immediately superficial to a palpable mammary tumor, and the collagen microenvironment was imaged in 8 and 12 weeks of age. Automated fiber extractions are shown overlaid on representative images from the 8- (a) and 12- (b) week time points. The bar graph (c) shows the ratio of the number of wavy fibers to total fibers found in the image. Fibers are labeled wavy if the distance along the fiber divided by the distance between fiber endpoints is greater than 1.08. Error bars indicate one standard deviation of the computed average wavy fractions among the nine images analyzed for each time point. Scale bar is 25 μm.

Thus, for perfectly straight fibers $W=1.0$ and for more wavy fibers $W>1.0$. We labeled a fiber as wavy if $W$ was greater than a threshold value of 1.08. Then, to compute the wavy fraction per image, the number of wavy fibers was divided by the total number of fibers found in each image. The resulting wavy-fraction values were averaged over all images at each time point and plotted in the bar graph shown in Fig. 8(c). We observe that the fraction of curvy fibers at the 8-week time point was approximately $0.19\pm 0.04$ and $0.75\pm 0.03$ at the 12-week time point. The error terms given here and error bars in Fig. 8(c) represent one standard deviation around the average wavy fractions for the nine images analyzed at each time point and indicate that there was close agreement between all images within a given time point.

4.

Discussion

In the present study, we compare preprocessing approaches prior to the application of the FIRE algorithm to identify fiber-level collagen characteristics in a series of SHG images of collagen in mammary tissue. Fiber extraction facilitates automated analysis of collagen features such as fiber number, length, and curvature. These features are important to researchers studying the role of the extracellular matrix in cancer progression. Computer-assisted interpretation of these fiber-level collagen patterns has the potential to generate more reliable and reproducible results compared to manual or transform-/filter-based quantification methods. Furthermore, an algorithm that identifies collagen fiber characteristics in tissue samples may enable large-scale studies of tumor-associated collagen signatures supporting the manual analysis performed previously.⁷

To our knowledge, FIRE has not been applied to SHG images of collagen in tissue. According to our testing, FIRE works well in some situations without any preprocessing or prefiltering. However, the algorithm fails when collagen fibers are densely packed or image quality is degraded, both of which are common occurrences while imaging collagen in tissue. Our work aims to extend FIRE’s application sphere to include complicated SHG images in tissue and to quantitatively compare the performance of a selection of preprocessing algorithms. Our results show that both the CT and the TF approaches are very promising and improve the fiber extraction accuracy achieved by the FIRE algorithm in many key situations. In addition, we demonstrate the application of our top-performing algorithm to extract in vivo fiber curvature changes during the development of a mouse mammary tumor. Although FIRE is used in our study for fiber extraction, other effective approaches that have been developed for vessel segmentation or neural diffusion mapping such as statistical tracking^40–42 may be effective in SHG image analysis. We believe the CT and TF methods would generally improve these algorithms as well.

Recent articles^43,44 have suggested that the CT method may be successfully applied in combination with other approaches for image processing such as fiber extraction, as we have demonstrated here. By selecting and thresholding the most representative scales, the CT-based method shows the best performance for both denoising the image and enhancing edge information, producing a better fiber extraction among the proposed preprocessing algorithms discussed in this article. Evidence supporting this claim is presented in Fig. 6, where the overall $F$-measure result was the highest for the CT method and notably a threefold improvement over a standard GF when considering all 25 test cases analyzed in this study. In addition, the CT-based method simplifies the often-difficult choice of selecting a threshold to binarize the image early in the FIRE process. Image thresholding can be difficult in low SNR and nonstationary images but may be alleviated through the application of more complicated thresholding techniques⁴⁵ or via the gray-level distance transform.⁴⁶ In our case, the inverse CT makes gray-level threshold selection simple by placing the background on the negative side of zero and the foreground on the positive side of zero, allowing the threshold to always remain at zero. Although gray-level thresholding is simplified, CT denoising adds the step of selecting how many curvelet coefficients to use. We chose a hard thresholding approach since it had robust performance for all cases we tested, though other soft thresholding or scale-adaptive thresholding techniques may be adopted to finely adjust the CT reconstruction. In addition, to take full advantage of the multiscale analysis of CT-based approaches, an optimal scale combination must be obtained according to the features of the images to take into account different fiber width, length, and dynamic intensity changes. The CurveAlign software¹⁶ previously developed in our group may be used to show the curvelet centers and directions of the fiber edges at a specified scale, which may be helpful for choosing the optimal scales and threshold of the curvelet coefficients.

The TF method produces slightly lower overall fiber segmentation accuracy compared to the CT method, as shown by the lower optimized $F$-measure score in Fig. 6. However, the optimized recall score of the TF method was higher than that of the CT method. This indicates the use of the TF method if recall is most critical; in other words, if priority is placed on not missing real fibers. In this study, we decided to balance recall and precision equally, therefore a missing fiber and a false alarm fiber were considered equally important (see $F$-measure calculation). When this is taken into account, although close, the CT method produced higher accuracy segmentations on average compared to the TF method.

The GF and SPTV methods produce similarly inferior segmentations compared to the CT and TF methods. One reason for this is that the GF and SPTV methods lack the ability to normalize the fibers in the images, such that bright or dim, or thick or thin fibers do not generate the same signal level in the output image. This lack of image normalization in the GF and SPTV methods causes difficulty in the threshold selection step. In addition, these filters do not enhance the ridges along the centers of the fibers, which is an attractive feature of both the CT and TF methods. The GF method is able to attenuate high-frequency noise but does not preserve edges. The SPTV method filters high-frequency noise and preserves edge information, however, it allows plateaus of high-signal level to remain in the image, such as those seen at the center of bright, thick fibers. For these reasons, the GF and SPTV methods are not ideal for preprocessing SHG images of collagen prior to fiber tracking.

Following our evaluation and comparison of these image-processing algorithms, we applied our top-performing algorithmic approach, termed CT-FIRE, to quantify collagen fiber shape changes over time in an in vivo mouse model for breast cancer. We observed a significantly larger fraction of highly curved fibers at the late time point compared to the early time point, indicating a quantifiable collagen matrix reorganization in the vicinity of the developing mammary tumor. Although the precise mechanisms underlying this observed matrix reorganization are currently unknown, we have demonstrated the power of a tool like CT-FIRE to quantify key aspects of these dynamic processes and others like it, enabling further studies of ECM remodeling.

This article focuses on image processing as a technique for quantifying structural information about collagen fibers in SHG images. However, it should be noted that there are a number of related techniques that use information about the polarization or directionality of the SHG signal to make inferences about collagen fiber orientation or estimates of the nonlinear susceptibility tensor.^31,47–52 These techniques have the potential drawback of adding costly components to the imaging system and often require multiple images to be captured per imaging frame, limiting their usefulness for applications such as in vivo imaging. Our goal here was to establish a robust technique for quantitative collagen architecture analysis of images captured with standard SHG imaging techniques.

It is worth mentioning, too, that although the CT and TF preprocessing methods can improve the results of the FIRE algorithm to some degree, they may do little about some intrinsic limitations of FIRE, such as the ability to properly segment crossing or cross-linked fibers, extremely curvy fibers, or fibers with gaps due to the fibers that travel in and out of the focal plane, as we observed in our testing. However, with the improvements provided by the combined approach of CT-FIRE, we anticipate being able to more accurately measure collagen fiber angle distributions in a highly automated fashion, thereby leading to better understanding of the interactions between the cells and collagen fibers. In order to link collagen architecture to cellular features, SHG imaging and CT-FIRE may be combined with complementary imaging techniques such as multiphoton-excited fluorescence imaging⁵³ and fluorescence lifetime imaging,⁵⁴ which allow imaging of both extrinsic and intrinsic fluorescences of tumor and stromal cells. In the future, accurate assessment of tumor–stromal interactions will help reveal prognosis or treatment response in diseases such as breast cancer.

5.

Conclusion

We demonstrate here an integrated approach for quantitative SHG collagen image analysis and algorithm evaluation. We show that the application of CT denoising as a preprocessing step for FIRE, a process we call CT-FIRE, performs more accurate fiber segmentations compared to other techniques we investigated in a variety of collagen images of human breast and mouse mammary tissues. We then demonstrate that CT-FIRE can automatically sense changes in collagen fiber curvature from images captured in an in vivo breast cancer mouse model. Our current work uses both MATLAB and Fiji³⁵ image-processing tools in combination. To make these approaches more widely accessible, we plan to develop a single Fiji plugin to perform the CT-FIRE process to produce 2-D and 3-D collagen fibers network extractions. Other future efforts will include the evaluation of multiple fiber-tracking algorithms applied to collagen fiber tracking in SHG images. Although this study focused solely on breast cancer, the use of these fiber quantification techniques should be easily adapted to SHG images of other collagen-related diseases. A MATLAB implementation of the CT-FIRE algorithm is available at http://loci.wisc.edu/software/ctfire.