breast cancer diagnosis and staging have been revolutionized by new molecular screening assays based on immunohistochemistry , fluorescence in situ hybridization , and reverse transcription polymerase chain reaction , which are all used to personalize care . however , there remains a significant group of breast cancer patients for whom these new techniques ultimately fail , due to several factors including varying patient genotype and primary or acquired resistance to drugs such as the her2/neu receptor targeting drug trastuzumab ( trade name herceptin ) . in addition , molecular screens are confounded by the high - degree of intratumor genetic diversity and often require extra tissue sections to be cut , stained and evaluated on top of the standard hematoxylin and eosin ( h&e ) preparation . new assays that predict patient outcome and response to treatment are therefore critically needed if we are to continue improving breast cancer treatment and prevention . one promising area of development is image based assays , which leverage high content imaging hardware and image analysis software to classify biological samples . in many cases , image based analysis does not require more than the standard histopathology h&e stained slides prepared as part of the normal clinical workflow . in this paper , we demonstrate the use of a new image - based assay for predicting patient outcome using information about tumor - stromal interactions from standard h&e stained histopathology specimens . the importance of stromal collagen in breast cancer is highlighted by the link between breast cancer , breast density , and the increased deposition of stromal collagen . interestingly , although mammographic density , which is attributable to collagen content , is one of the largest risk factors for the development of breast tumors , there is currently no clinical intervention based on mammographic density alone . this is due in part to the lack of a clear correlation observed between increased mammographic density and patient outcome . most of the work to date has defined mammographic density as a etiological factor and not as a prognostic factor . they reported that women with intermediate and high breast density had a significantly elevated risk to develop a local breast cancer recurrence . however , follow - up clinical trials that incorporate additional risk factors such as obesity are needed to examine the possible prognostic value of mammographic density in large and diverse patient cohorts before using density as a possible clinical target . as well recently there has been an effort to investigate the underlying contributor to mammographic density , focusing on one of the largest components present in the dense stroma , collagen . several studies have shown a link between collagen remodeling and the invasion and progression of mammary cancer in mouse models . furthermore , there was a link observed between collagen morphology , particularly collagen alignment , and breast cancer patient outcome . first introduced the so called tumor associated collagen signature ( tacs ) nomenclature to describe collagen alignment patterns . tacs-2 is observed as straightened fibers aligned tangentially around developing tumors , while tac-3 is seen as radially aligned fibers that facilitate local invasion . conklin et al . qualitatively searched for these patterns in human breast cancer samples and through extensive manual analysis found that the presence of the tacs-3 alignment phenotype was a prognostic indicator for disease free and disease specific survival ( dfs and dss respectively ) for invasive breast cancer patients . our quantitative study , presented here , computationally builds on this previous work by defining an algorithmic model for tacs-3 and applying this model to the same cohort of patients . previous collagen alignment studies have largely been facilitated by the development of second harmonic generation ( shg ) microscopy techniques , which have the ability to capture high contrast images of the collagen fiber extracellular matrix without the need for exogenous stains . for example , changes in the ratio of the forward shg ( fshg ) to backward propagating shg signal have been recently linked to breast tumor progression and positive lymph node status . shg directionality was also used by ajeti et al . to quantify the collagen composition in breast cancer models , while ambekar et al . used fourier transform and polarization - resolved shg imaging to differentiate malignant from benign tissues in breast biopsies . in addition , many new computational techniques are being developed to quantify patterns observed in shg images . provides three - dimensional orientation and radius information about fibers in shg images . due to the fibrous nature of the collagen matrix , shg images are particularly well - suited for the curvelet transform ( ct ) , which is a multiscale , orientation sensitive version of the wavelet transform . the ct and combined fiber tracking methods have been applied to extract fiber orientation , length , curvature and radius from shg images of collagen . one key feature however that is missing from all of the available image analysis techniques is the ability to incorporate cellular information into the analysis . the interaction between tumor cells and collagen fibers can not be fully assessed without integration of information about cellular morphology and associated collagen morphology . as well this information is critical for finding regions of interest ( roi ) with tacs , an essential task for any type of high - throughput screening where manual searching is not practical . herein , we describe a computational protocol that achieves this goal by integrating information about collagen fibers from shg images with information about cells captured through bright field imaging of standard h&e stained slides to perform highly automated , prognostic tacs-3 scoring . in order for tacs to become a useful and fully validated biomarker , it must be screened for in several large studies containing many patients and diverse populations . in addition , besides screening in heterogeneous populations , it ideally needs to be screened in diverse sample types to account for possible subtle differences in surgery , pathology or sample preparation that could negatively impact sample consistency . this ability to rapidly screen in many sample types of large diverse populations would also open the door for tacs to be explored in other cancer types such as pancreatic and renal cancer . heretofore , there has not been a method that automates enough of the process to enable such large scale adaptation . in previous studies , collagen fiber angles were measured by hand , one at a time , using imagej roi marking tools . these experiments used information gathered a priori or from autofluorescence to identify tumor - stromal boundaries . conklin et al . manually captured each individual image , used bright field images to manually identify tumor - stromal boundaries , and manually estimated collagen fiber angles . in each of these cases , many subjective decisions were made while identifying which areas to image , which fibers to measure and what should be considered a tumor - stromal boundary . there has been progress made in automating the fiber angle analysis steps of this task . however , none of these methods can automate all four steps of the tacs analysis process , which are : ( 1 ) image capture , ( 2 ) fiber angle measurement , ( 3 ) tumor - stromal boundary identification , and ( 4 ) relative angle measurement between fiber and boundary . in this paper , we use image analysis and supervised learning techniques to enable the automation of each of these tasks . the block diagram of our imaging and analysis protocol is shown in figure 1 . starting with the previously - imaged invasive breast cancer tissue microarray ( tma ) , we captured registered , whole - slide shg and bright field images , extracted fibers from the shg images , identified tumor - stromal boundaries from the bright field images , and measured relative angles , all in a scripted pipelined process that requires little human intervention . we believe that this method will allow significantly larger scale studies to be performed in order to validate tacs-3 as a prognostic biomarker in breast cancer and potentially other cancer types , and to investigate if tacs-3 can be used to predict patient response to targeted therapies . the red steps are performed with wiscscan , the blue steps are performed with imagej / fiji , and the green steps are performed with matlab based tools . the left side shows the steps performed on the second harmonic generation images , while the right side shows the steps performed on the bright field images . the tma used here was the same as that used by conklin et al . for the manual collagen alignment analysis . the clinical profiles of all patients whose tissue was included in this tma have been described in a previous study . all tissue and patient information used in this study were acquired following institutional review board approval . tumor tissues from 353 patients diagnosed with invasive carcinoma were resected by the same surgeon between 1981 and 1995 . after tumors smaller than 5 mm and severely damaged samples were excluded , 196 patients remained for analysis . sections of 4 m thickness were cut from archived tma blocks containing 1.0 mm diameter tissue cores , placed on glass slides , stained with h&e and mounted under a glass coverslip . patients were followed for a median of 6.2 years , ranging from 1 month to 18.6 years . all samples in this study were imaged with the custom built integrated fshg / bright field imaging system shown in figure 2 . a mira 900 ti : sapphire laser ( coherent , santa clara , ca ) tuned to 780 nm , with a pulse length of approximately 100 fs , was directed through a pockel 's cell ( conoptics , danbury , ct , usa ) , half and quarter waveplates ( thorlabs , newton , nj , usa ) , beam expander ( thorlabs ) , a 3 mm galvanometer driven mirror pair ( cambridge , bedford , ma ) , a scan / tube lens pair ( thorlabs ) , through a dichroic beam splitter ( semrock , rochester , ny ) and focused by a 20x/0.75na objective ( nikon , melville , ny ) . shg light was collected in the forward direction with a 0.54 na condenser ( thorlabs ) and filtered with an interference filter centered at 390 nm with a full width at half maximum bandwidth of 22.4 nm ( semrock ) . the back aperture of the condenser lens was imaged onto the 5 mm aperture of a 7422 - 40p photomultiplier tube ( hamamatsu , hamamatsu , japan ) the signal from which was amplified with a c7319 integrating amplifier ( hamamatsu ) and sampled with an analog to digital converter ( innovative integration , simi valley , ca ) . timing between the galvo scanners , signal acquisition , and motorized stage positioning was achieved using our custom software called wiscscan . the rapid automated modular microscope system ( applied scientific instrumentation , eugene , or ) served as our microscope base and we used asi motorized translation stages for x , y , and z motion control . the shg light source was verified to be circularly polarized at the sample using the protocol of chen et al . shg images were captured as stacks of three images spaced 3 m apart , then z - projected to improve field flatness . bright field images were captured with the same system using a mcwhl2 white led lamp ( thorlabs ) set up for kohler illumination . white light from this lamp was separated from shg light traveling through the condenser assembly using a short pass dichroic mirror with a cutoff at 670 nm ( semrock ) . an rgb camera ( qimaging , surrey , bc , canada ) was used to capture bright field images through wiscscan to allow for acquisition within a single application . both shg and white light images were tiled with 10% overlap using automation provided by wiscscan . stage positions for individual images and pixel size data were stored in bio - formats image metadata and this was then used by the grid / collection stitching imagej plugin to reassemble a high - resolution large field of view image of the entire tma . when capturing large field of view images , the sample plane often walks out of the in - focus imaging plane as the stage is translated over large distances in x or y. we alleviated this issue using the continuous reflection interface sampling and positioning autofocus system ( applied scientific instrumentation ) , which maintained an accurate distance between the coverslip and the objective throughout the whole slide stitched image capture . this allowed for a single bright field image to be captured at each location rather than a z - stack , improving capture speed , reconstruction speed , and reducing production of unnecessary data . after shg and white light images were captured and stitched , the two modalities were manually registered with the landmark correspondences imagej plugin using five control points per image . the image of the entire tma was registered in a single step , and then each individual tma core was cropped out of the full tma image , producing 196 images . the resulting tma core images were each 2048  2048 pixels , consisting of four eight - bit channels . the first three channels represented the red , green and blue planes of the white light image , while the fourth channel contained the shg information . photo of imaging system ( a ) and optical block diagram panel ( b ) . pmt = photomultiplier tube , qcam = qimaging rgb camera , led = bright field lamp , dbs = dichroic beamsplitter , tl = tube lens , cl = condenser lens , l = lens , m = mirror , be = beam expander , z = z - direction translation , xy = xy translation , xyg = xv galvonometer driven mirrors , pc = pockel 's cell , ramm = rapid automated modular microscope ( asi ) , 2 , 4 = half and quarter waveplates our tacs-3 model was based on previously published observations relating collagen structure to breast cancer progression and survival . in these studies , the first step in the tacs-3 scoring process was the identification of groups of straightened , aligned collagen fibers . the second step was to determine if those fibers terminate at or near regions of epithelial cells at steep angles . if a fiber met both of these criteria , then it was considered tacs-3 positive . if one or more tacs-3 positive fibers were found in a sub - region of an image , then that region was scored tacs-3 positive . the number of regions with tacs-3 positive scores was then used to score the entire image . there were many details in these steps and defining parameters to account for each step would have produced a potentially fragile model . instead , we have implemented a supervised learning approach that allows the data to most appropriately define the model . the first classifier was trained to find epithelial regions in the images using a small training set of annotated rois . the resultant epithelial cell model was then used to segment epithelial cell regions within the entire cohort of images . features describing the epithelial regions were then combined with features derived from our fiber extraction algorithm and were fed into a second classifier , which was trained to score each image as being tacs-3 positive or negative based on a training set of annotated images . tacs-3 scores were then fed into a cox proportional hazard model to regress to censored survival data . our entire image cohort can be represented by a set s of registered shg and bright field ( rgb ) images . each image ii( , g)s was composed of pixels =(xi , yi)r,0<xi , yi < n and the function g( ) which mapped each pixel to a quartet of intensity values corresponding to the r , g , b , and shg intensity channels of the image . the first step of our protocol involved the extraction of collagen fiber objects from the shg channel . next , epithelial cell clusters were segmented from the r , g , and b channels . each fiber from the first step was then associated with epithelial cell clusters from the second step to create a feature set for each fiber . average feature values for all fibers in the image were then used to classify images based on a training set of images . we applied a technique called ct - fire to the shg images to enhance , trace and extract a network of collagen fibers for each shg image i( , gshg ) . ct - fire combines the advantage of the ct for denoising the image and enhancing the fiber ridge features with the advantage of a fiber tracing algorithm for automatic fiber extraction , being capable of extracting fiber geometric information such as length , angle , width , and curvature of each fiber . we applied the fast discrete ct ( fdct ) to capture a collection of coefficients c in curvelet space , which are defined as the inner product of the input shg image channel with each of the curvelet basis functions . where   is the digital curvelet waveform and jkl represent the scale , orientation , and location indices , respectively . we used the open source fdct matlab ( the mathworks , natick , ma , usa ) library and specifically the  wrapping  version of the fdct due to its simplicity . to denoise the image , we set all curvelet coefficients to zero that fall below a user defined threshold t as shown below this threshold was determined empirically on a small subset of shg images to determine the appropriate level of noise reduction . the inverse fdct was then applied to reconstruct an edge enhanced , noise reduced version of the shg image . after reconstruction , ct - fire traced fibers , using the method of stein et al . the distance transform computed the distance from each foreground pixel to the nearest background pixel . fiber branches were formed by creating regions surrounding each local maxima , the size of which were defined by the result of the distance transform at the location of the local maximum point . this process was repeated until no new local maxima were found indicating the end of a fiber branch . short branches were then pruned from the network and closely spaced , similarly oriented fibers were merged . fiber width ( fw ) was quantified for each extracted fiber by averaging the fiber widths ( 2ri ) at n points that were used to form the fiber where ri is the fiber radius at the i point , estimated by the result of the distance transform at that location . fiber straightness ( fs ) was quantified for each extracted fiber by dividing the distance between the end points of the fiber ( dn ) by the distance along the path of the fiber ( d0 ) . thus for perfectly straight fibers fs=1.0 and wavy fibers fs<1.0 . after fiber objects have been extracted from each of the images , we next segment epithelial cell regions . the tacs-3 phenotype consists of straightened aligned collagen fibers that terminate near regions of epithelial cells such that the angles of the collagen fibers appear perpendicular to the epithelial stromal boundary . detecting this tacs-3 phenotype requires knowledge of the locations of epithelial cells within the sample . we must then identify regions of epithelial cell clusters and identify a boundary between the epithelial cells and surrounding stroma . step 1 used the trainable weka segmentation imagej plugin to find epithelial cell nuclei and step 2 applied a cascaded matched filter , threshold operation to identify clusters and boundaries . for step 1 , a training set of 15 cropped 256  256 pixel images denoted as ti(t , g)s was created that contains representative features from five classes : epithelial cell nuclei , other cell nuclei ( including lymphocytes and fibroblasts ) , cytoplasm , collagen , and background . a further subgroup of pixels within the training images at were annotated as belonging to each class wk(a),k(1,  ,5 ) . a feature vector pi(t , g ) was computed for each pixel and each channel of the training image where i ( l  .d ) was the feature index and d was the dimensionality of the feature subspace . the feature set we used is listed in table 1 and incorporates features at five scales for a total of 80 feature planes for each image channel . the detailed implementations for each of these features are given in the online documentation for the weka segmentation plugin . we then used a multithreaded implementation of the random forests classifier with a forest of 200 trees and two random features per node to build a model based on a and wk . the trained model was then applied to every pixel  in the cohort producing a probability map for each class and each image using a scripted version of the plugin . features used in epithelial cell segmentation and tacs-3 fiber classification tasks for the second step in the segmentation process , the epithelial class probability map was filtered with a gaussian filter matched to the average width of the epithelial cell nuclei ( three microns ) and thresholded such that the top 80% of resulting pixels were retained . the resulting image was then filtered with a gaussian filter matched to the width of the average sized epithelial cell cluster ( 25 microns ) , then finally thresholded such that , again , the top 80% of resulting pixels were retained . following the final threshold step , regions smaller than 50 pixels in area were discarded and a mask was generated with epithelial cell clusters in the foreground and all else in the background . epithelial mask pixels are represented here as ei while epithelial region boundary pixels were created using an eight - connected neighborhood and are denoted as bi . mask images were saved as tiff files and read , along with the extracted fiber data , into the custom , open source curvealign software , described more below , for fiber / epithelial region feature extraction . outlines of the resulting mask files were overlaid onto the original white light images to qualitatively validate the segmentation accuracy of the applied epithelial region model . in the sections above , we described our methods for epithelial cluster segmentation and collagen fiber extraction . with these two pieces of information , we associated fibers with epithelial cell clusters and measured the interaction between the two using the features described here . this tool started by reading in a fiber database file ( generated by ct - fire ) and an epithelial mask file ( generated by our epithelial segmentation script ) . a feature vector pi was then built for each fiber endpoint vir in the image . the feature vector was populated directly with features derived above in the fiber extraction section including fiber length , curvature , radius and grey level . many of the features in this section rely heavily on the nearest neighbor search routine which is formulated here as where d={x1,  ,xn } is a set of vectors in r , q is a query vector , (x , q ) is the euclidean norm x  q2 , and the result (x , q ) is a vector of n points in x that are nearest to each point in q. given a collection of points on a two - dimensional plane , if we select a query point , this algorithm will return the nearest neighbors within our collection of points . for example , a fiber end point can be used as the query point and we can use this algorithm to search for the nearest point in the list of epithelial cell boundary points , as described below . in addition , some of the features compute a metric for alignment using vector addition according to the following algorithm where {0,  , } is a vector of orientations associated with the vector of n positions in x. the factor of two is included since we used fiber orientations supported from 0to rather than full 0to2 oriented direction vectors . in words , the alignment metric is calculated as the normalized vector sum of orientation vectors . the larger the vector sum , the more aligned the group of fibers . on the other hand , if the vector sum is small , then the group of fibers is more randomly oriented . fiber density features were computed as the average distance from the current fiber endpoint to the n = 2 , 4 , 8 , and 16 nearest neighbors . the density features for fiber endpoint v are therefore given by : if fiber density is higher , then this result will be lower , since it is measuring the average distance between fibers . fiber alignment features were computed as the absolute values of the vector sum of the nearest neighbor fiber endpoints and are given by   for fiber endpoint . in the equation for fan , v,[. . ] features that incorporate epithelial cell information included distance to nearest epithelial mask point dev = (v(ei , v ) ) and distance to nearest epithelial region boundary dbv = (v(biv ) ) . these features had the same value if the fiber end point was outside an epithelial cell region ; however , they were different if the end point was colocal with an epithelial region . angle with respect to nearest epithelial region boundary point was computed as abv=([v,(bi , v ] ) and angle with respect to nearest  extension boundary intersection point  was given by , where the set of points in l was computed by taking the intersection of all epithelial boundary points bi and a line of length 2r extending from the fiber endpoint at an angle  and is formulated below where bres indicates a modified bresenham algorithm which is used to find all pixels along a line between two points . the term q=[rexp(v ) ] is the offset from v in the x and y directions . for this last feature , three values of r ( 50 , 100 and 200 m ) were calculated . these three lengths corresponded to 5 , 10 , and 20 times the diameter of a typical epithelial cell and were selected based on estimates of intercellular signaling distances . if no intersection was found , then the   feature value was set to zero . the angle of the tumor stromal boundary line   was estimated by fitting a quadratic to nine contiguous points on the boundary surrounding the intersection point ( or nearest boundary point in the previous feature ) and computing the tangent angle of the line fit at the midpoint . the steps in the process of relative angle feature extraction are diagrammed in figure 3 . panel a shows a single fiber ( green ) and epithelial region boundary ( red ) with one highlighted fiber endpoint . zoomed versions of panel a are shown in panels b , c and d where individual image pixels are represented as filled circles . the nearest distance from to the boundary is indicated in panel b , intersection between the endpoint extension line and the boundary is shown in panel c and the quadratic curve fit to the boundary at the intersection point and tangent line are shown in panel d each of these fiber level features pi were calculated for every fiber endpoint vi in the cohort . fiber level features were then averaged among all fibers in a given image and training was performed with a subset of 16 images itii that had been manually annotated as being tacs-3 positive or negative . a linear support vector machine ( svm ) was used to build a model , which was then applied to all 196 images in the cohort for classifying each image as being tacs-3 positive or negative . the tacs-3 scores were correlated with dfs and dss data using the cox - proportional - hazards regression method . dfs was defined as the time from date of diagnosis to the first date of recurrence and dss was defined as the time from diagnosis to death from breast cancer or date of last follow - up evaluation . in both cases , all other events were censored . the kaplan - meier method was used to compare dfs and dss between tacs-3 negative and tacs-3 positive patients . correlations between manual and computationally generated tacs-3 scores were made using the pearson 's linear correlation coefficient . the tma used here was the same as that used by conklin et al . for the manual collagen alignment analysis . the clinical profiles of all patients whose tissue was included in this tma have been described in a previous study . all tissue and patient information used in this study were acquired following institutional review board approval . tumor tissues from 353 patients diagnosed with invasive carcinoma were resected by the same surgeon between 1981 and 1995 . after tumors smaller than 5 mm and severely damaged samples were excluded , 196 patients remained for analysis . sections of 4 m thickness were cut from archived tma blocks containing 1.0 mm diameter tissue cores , placed on glass slides , stained with h&e and mounted under a glass coverslip . patients were followed for a median of 6.2 years , ranging from 1 month to 18.6 years . all samples in this study were imaged with the custom built integrated fshg / bright field imaging system shown in figure 2 . a mira 900 ti : sapphire laser ( coherent , santa clara , ca ) tuned to 780 nm , with a pulse length of approximately 100 fs , was directed through a pockel 's cell ( conoptics , danbury , ct , usa ) , half and quarter waveplates ( thorlabs , newton , nj , usa ) , beam expander ( thorlabs ) , a 3 mm galvanometer driven mirror pair ( cambridge , bedford , ma ) , a scan / tube lens pair ( thorlabs ) , through a dichroic beam splitter ( semrock , rochester , ny ) and focused by a 20x/0.75na objective ( nikon , melville , ny ) . shg light was collected in the forward direction with a 0.54 na condenser ( thorlabs ) and filtered with an interference filter centered at 390 nm with a full width at half maximum bandwidth of 22.4 nm ( semrock ) . the back aperture of the condenser lens was imaged onto the 5 mm aperture of a 7422 - 40p photomultiplier tube ( hamamatsu , hamamatsu , japan ) the signal from which was amplified with a c7319 integrating amplifier ( hamamatsu ) and sampled with an analog to digital converter ( innovative integration , simi valley , ca ) . timing between the galvo scanners , signal acquisition , and motorized stage positioning was achieved using our custom software called wiscscan . the rapid automated modular microscope system ( applied scientific instrumentation , eugene , or ) served as our microscope base and we used asi motorized translation stages for x , y , and z motion control . the shg light source was verified to be circularly polarized at the sample using the protocol of chen et al . shg images were captured as stacks of three images spaced 3 m apart , then z - projected to improve field flatness . bright field images were captured with the same system using a mcwhl2 white led lamp ( thorlabs ) set up for kohler illumination . white light from this lamp was separated from shg light traveling through the condenser assembly using a short pass dichroic mirror with a cutoff at 670 nm ( semrock ) . an rgb camera ( qimaging , surrey , bc , canada ) was used to capture bright field images through wiscscan to allow for acquisition within a single application . both shg and white light images were tiled with 10% overlap using automation provided by wiscscan . stage positions for individual images and pixel size data were stored in bio - formats image metadata and this was then used by the grid / collection stitching imagej plugin to reassemble a high - resolution large field of view image of the entire tma . when capturing large field of view images , the sample plane often walks out of the in - focus imaging plane as the stage is translated over large distances in x or y. we alleviated this issue using the continuous reflection interface sampling and positioning autofocus system ( applied scientific instrumentation ) , which maintained an accurate distance between the coverslip and the objective throughout the whole slide stitched image capture . this allowed for a single bright field image to be captured at each location rather than a z - stack , improving capture speed , reconstruction speed , and reducing production of unnecessary data . after shg and white light images were captured and stitched , the two modalities were manually registered with the landmark correspondences imagej plugin using five control points per image . the image of the entire tma was registered in a single step , and then each individual tma core was cropped out of the full tma image , producing 196 images . the resulting tma core images were each 2048  2048 pixels , consisting of four eight - bit channels . the first three channels represented the red , green and blue planes of the white light image , while the fourth channel contained the shg information . photo of imaging system ( a ) and optical block diagram panel ( b ) . pmt = photomultiplier tube , qcam = qimaging rgb camera , led = bright field lamp , dbs = dichroic beamsplitter , tl = tube lens , cl = condenser lens , l = lens , m = mirror , be = beam expander , z = z - direction translation , xy = xy translation , xyg = xv galvonometer driven mirrors , pc = pockel 's cell , ramm = rapid automated modular microscope ( asi ) , 2 , 4 = half and quarter waveplates our tacs-3 model was based on previously published observations relating collagen structure to breast cancer progression and survival . in these studies , the first step in the tacs-3 scoring process was the identification of groups of straightened , aligned collagen fibers . the second step was to determine if those fibers terminate at or near regions of epithelial cells at steep angles . if a fiber met both of these criteria , then it was considered tacs-3 positive . if one or more tacs-3 positive fibers were found in a sub - region of an image , then that region was scored tacs-3 positive . the number of regions with tacs-3 positive scores was then used to score the entire image . there were many details in these steps and defining parameters to account for each step would have produced a potentially fragile model . instead , we have implemented a supervised learning approach that allows the data to most appropriately define the model . the first classifier was trained to find epithelial regions in the images using a small training set of annotated rois . the resultant epithelial cell model was then used to segment epithelial cell regions within the entire cohort of images . features describing the epithelial regions were then combined with features derived from our fiber extraction algorithm and were fed into a second classifier , which was trained to score each image as being tacs-3 positive or negative based on a training set of annotated images . tacs-3 scores were then fed into a cox proportional hazard model to regress to censored survival data . our entire image cohort can be represented by a set s of registered shg and bright field ( rgb ) images . each image ii( , g)s was composed of pixels =(xi , yi)r,0<xi , yi < n and the function g( ) which mapped each pixel to a quartet of intensity values corresponding to the r , g , b , and shg intensity channels of the image . the first step of our protocol involved the extraction of collagen fiber objects from the shg channel . next , epithelial cell clusters were segmented from the r , g , and b channels . each fiber from the first step was then associated with epithelial cell clusters from the second step to create a feature set for each fiber . average feature values for all fibers in the image were then used to classify images based on a training set of images . we applied a technique called ct - fire to the shg images to enhance , trace and extract a network of collagen fibers for each shg image i( , gshg ) . ct - fire combines the advantage of the ct for denoising the image and enhancing the fiber ridge features with the advantage of a fiber tracing algorithm for automatic fiber extraction , being capable of extracting fiber geometric information such as length , angle , width , and curvature of each fiber . we applied the fast discrete ct ( fdct ) to capture a collection of coefficients c in curvelet space , which are defined as the inner product of the input shg image channel with each of the curvelet basis functions . where   is the digital curvelet waveform and jkl represent the scale , orientation , and location indices , respectively . we used the open source fdct matlab ( the mathworks , natick , ma , usa ) library and specifically the  wrapping  version of the fdct due to its simplicity . to denoise the image , we set all curvelet coefficients to zero that fall below a user defined threshold t as shown below this threshold was determined empirically on a small subset of shg images to determine the appropriate level of noise reduction . the inverse fdct was then applied to reconstruct an edge enhanced , noise reduced version of the shg image . after reconstruction , ct - fire traced fibers , using the method of stein et al . the distance transform computed the distance from each foreground pixel to the nearest background pixel . fiber branches were formed by creating regions surrounding each local maxima , the size of which were defined by the result of the distance transform at the location of the local maximum point . this process was repeated until no new local maxima were found indicating the end of a fiber branch . short branches were then pruned from the network and closely spaced , similarly oriented fibers were merged . fiber width ( fw ) was quantified for each extracted fiber by averaging the fiber widths ( 2ri ) at n points that were used to form the fiber where ri is the fiber radius at the i point , estimated by the result of the distance transform at that location . fiber straightness ( fs ) was quantified for each extracted fiber by dividing the distance between the end points of the fiber ( dn ) by the distance along the path of the fiber ( d0 ) . after fiber objects have been extracted from each of the images , we next segment epithelial cell regions . the tacs-3 phenotype consists of straightened aligned collagen fibers that terminate near regions of epithelial cells such that the angles of the collagen fibers appear perpendicular to the epithelial stromal boundary . detecting this tacs-3 phenotype requires knowledge of the locations of epithelial cells within the sample . we must then identify regions of epithelial cell clusters and identify a boundary between the epithelial cells and surrounding stroma . step 1 used the trainable weka segmentation imagej plugin to find epithelial cell nuclei and step 2 applied a cascaded matched filter , threshold operation to identify clusters and boundaries . for step 1 , a training set of 15 cropped 256  256 pixel images denoted as ti(t , g)s was created that contains representative features from five classes : epithelial cell nuclei , other cell nuclei ( including lymphocytes and fibroblasts ) , cytoplasm , collagen , and background . a further subgroup of pixels within the training images at were annotated as belonging to each class wk(a),k(1,  ,5 ) . a feature vector pi(t , g ) was computed for each pixel and each channel of the training image where i ( l  .d ) was the feature index and d was the dimensionality of the feature subspace . the feature set we used is listed in table 1 and incorporates features at five scales for a total of 80 feature planes for each image channel . the detailed implementations for each of these features are given in the online documentation for the weka segmentation plugin . we then used a multithreaded implementation of the random forests classifier with a forest of 200 trees and two random features per node to build a model based on a and wk . the trained model was then applied to every pixel  in the cohort producing a probability map for each class and each image using a scripted version of the plugin . features used in epithelial cell segmentation and tacs-3 fiber classification tasks for the second step in the segmentation process , the epithelial class probability map was filtered with a gaussian filter matched to the average width of the epithelial cell nuclei ( three microns ) and thresholded such that the top 80% of resulting pixels were retained . the resulting image was then filtered with a gaussian filter matched to the width of the average sized epithelial cell cluster ( 25 microns ) , then finally thresholded such that , again , the top 80% of resulting pixels were retained . following the final threshold step , regions smaller than 50 pixels in area were discarded and a mask was generated with epithelial cell clusters in the foreground and all else in the background . epithelial mask pixels are represented here as ei while epithelial region boundary pixels were created using an eight - connected neighborhood and are denoted as bi . mask images were saved as tiff files and read , along with the extracted fiber data , into the custom , open source curvealign software , described more below , for fiber / epithelial region feature extraction . outlines of the resulting mask files were overlaid onto the original white light images to qualitatively validate the segmentation accuracy of the applied epithelial region model . in the sections above , we described our methods for epithelial cluster segmentation and collagen fiber extraction . with these two pieces of information , we associated fibers with epithelial cell clusters and measured the interaction between the two using the features described here . this tool started by reading in a fiber database file ( generated by ct - fire ) and an epithelial mask file ( generated by our epithelial segmentation script ) . a feature vector pi was then built for each fiber endpoint vir in the image . the feature vector was populated directly with features derived above in the fiber extraction section including fiber length , curvature , radius and grey level . many of the features in this section rely heavily on the nearest neighbor search routine which is formulated here as where d={x1,  ,xn } is a set of vectors in r , q is a query vector , (x , q ) is the euclidean norm x  q2 , and the result (x , q ) is a vector of n points in x that are nearest to each point in q. given a collection of points on a two - dimensional plane , if we select a query point , this algorithm will return the nearest neighbors within our collection of points . for example , a fiber end point can be used as the query point and we can use this algorithm to search for the nearest point in the list of epithelial cell boundary points , as described below . in addition , some of the features compute a metric for alignment using vector addition according to the following algorithm where {0,  , } is a vector of orientations associated with the vector of n positions in x. the factor of two is included since we used fiber orientations supported from 0to rather than full 0to2 oriented direction vectors . in words , the larger the vector sum , the more aligned the group of fibers . on the other hand , if the vector sum is small , then the group of fibers is more randomly oriented . fiber density features were computed as the average distance from the current fiber endpoint to the n = 2 , 4 , 8 , and 16 nearest neighbors . the density features for fiber endpoint v are therefore given by : if fiber density is higher , then this result will be lower , since it is measuring the average distance between fibers . fiber alignment features were computed as the absolute values of the vector sum of the nearest neighbor fiber endpoints and are given by   for fiber endpoint . in the equation for fan , v,[. . ] features that incorporate epithelial cell information included distance to nearest epithelial mask point dev = (v(ei , v ) ) and distance to nearest epithelial region boundary dbv = (v(biv ) ) . these features had the same value if the fiber end point was outside an epithelial cell region ; however , they were different if the end point was colocal with an epithelial region . angle with respect to nearest epithelial region boundary point was computed as abv=([v,(bi , v ] ) and angle with respect to nearest  extension boundary intersection point  was given by , where the set of points in l was computed by taking the intersection of all epithelial boundary points bi and a line of length 2r extending from the fiber endpoint at an angle  and is formulated below where bres indicates a modified bresenham algorithm which is used to find all pixels along a line between two points . the term q=[rexp(v ) ] is the offset from v in the x and y directions . for this last feature , three values of r ( 50 , 100 and 200 m ) these three lengths corresponded to 5 , 10 , and 20 times the diameter of a typical epithelial cell and were selected based on estimates of intercellular signaling distances . if no intersection was found , then the   feature value was set to zero . the angle of the tumor stromal boundary line   was estimated by fitting a quadratic to nine contiguous points on the boundary surrounding the intersection point ( or nearest boundary point in the previous feature ) and computing the tangent angle of the line fit at the midpoint . panel a shows a single fiber ( green ) and epithelial region boundary ( red ) with one highlighted fiber endpoint . zoomed versions of panel a are shown in panels b , c and d where individual image pixels are represented as filled circles . the nearest distance from to the boundary is indicated in panel b , intersection between the endpoint extension line and the boundary is shown in panel c and the quadratic curve fit to the boundary at the intersection point and tangent line are shown in panel d each of these fiber level features pi were calculated for every fiber endpoint vi in the cohort . fiber level features were then averaged among all fibers in a given image and training was performed with a subset of 16 images itii that had been manually annotated as being tacs-3 positive or negative . a linear support vector machine ( svm ) was used to build a model , which was then applied to all 196 images in the cohort for classifying each image as being tacs-3 positive or negative . the tacs-3 scores were correlated with dfs and dss data using the cox - proportional - hazards regression method . dfs was defined as the time from date of diagnosis to the first date of recurrence and dss was defined as the time from diagnosis to death from breast cancer or date of last follow - up evaluation . in both cases , the kaplan - meier method was used to compare dfs and dss between tacs-3 negative and tacs-3 positive patients . correlations between manual and computationally generated tacs-3 scores were made using the pearson 's linear correlation coefficient . registered shg and bright field images of a subsample of the tma are shown in figure 4 along with two zoomed versions of regions within the image . shg information is added as an alpha channel on top of the raw rgb bright field image and pseudo colored yellow . the fully zoomed panel shows the detail available in the full resolution images captured with the 20x , 0.75 na lens and shows a region with a positive tacs-3 signature . a collection of three more tacs-3 positive and three tacs-3 negative regions were cropped out of the tma images and shown in figure 5 . these images illustrate the features that are common to the tacs-3 signature including straightened , aligned fibers terminating in or near regions of epithelial cells at near perpendicular angles with respect to the epithelial region border . in addition , the tacs-3 negative cases show wavy fibers , fibers that terminate at adipose tissue , and a curved fiber encapsulating an epithelial cell cluster [ figure 5d  f , respectively ] . the forward second harmonic generation image was overlaid upon the bright field image for an entire h&e stained tissue microarray slide . scale bar = 25 m examples of tumor associated collagen signatures ( tacs-3 ) positive and negative regions within our cohort of breast cancer samples . collagen fibers ( yellow ) registered and overlaid onto bright field images of an h&e stained tissues . these images illustrate the features of the tacs-3 signature ( green arrows , top row ) , particularly , straightened aligned collagen fibers that terminate at steep angles relative to epithelial region boundaries ( a - c ) . three tacs-3 negative cases are shown for comparison ( red arrows , bottom row ) , where collagen fibers are wavy ( d ) , terminate at adipose tissue ( e ) , and wrap around epithelium ( f ) scale bar = 25 m a sample of our fiber extraction and epithelial region segmentation results are shown in figures 6 and 7 , respectively . in both cases , epithelial region segmentation and fiber extraction were observed to accurately represent the data . the orientations of the epithelial cell region boundaries were compared to collagen fiber angles derived from the results of a fiber object extraction algorithm ct - fire which has been shown to perform well in comparison to other techniques . a representative sample of the results produced by this algorithm are shown in figure 6 . the intermediate product after the ct denoising step is shown in figure 6a , while the extracted fiber network is shown overlaid on the original shg image as shown in figure 6b . although some fibers are over - or under - segmented ( annotated by green arrows ) , most of the extracted fibers properly represent the data . figure 7 clearly demonstrates the ability of our epithelial cell segmentation algorithm to properly classify many of the regions of epithelial cells as positive . however , a few small regions of stromal fibroblasts and endothelial cells are included in the epithelial cell regions ( annotated by green arrows ) . although these errors occurred occasionally throughout the cohort , the noise they generated did not overcome the tacs-3 signal . the boundary smoothness was dependent on the selection of our filter widths and binary mask thresholds . these parameters were selected to accurately represent the boundary orientation at the spatial scale of the epithelial cell regions . sample fiber extraction results . the resulting image after curvelet denoising ( a ) shows likely fiber pixels in white and likely background pixels in grey . the extracted fiber network is overlaid on the original second harmonic generation image ( b ) with many appropriate segmentations and a few under - and over - segmentations ( green arrows ) . the raw probability map produced by the trainable weka segmentation imagej plugin ( a ) is filtered to estimate epithelial cluster density ( b ) and thresholded ( c ) to produce epithelial region boundaries which are overlaid onto original bright field images to validate the segmentation ( d ) scale bar = 100 m although correlation with survival is our ultimate goal , automated tacs-3 scores should also correlate with manual scores for each of the images . the pearson linear correlation coefficient was used to determine this correlation , the results of which are tabulated in table 2 . score 1 was the number of tacs-3 positive regions divided by the total number of regions analyzed , score 2 was the average number of tacs-3 positive votes per region among three observers , and score 3 indicated if one or more region received a tacs-3 positive rating . table 2 shows positive correlation between all manual scoring methods and our computational scoring system presented here , with the highest correlation observed to be with manual score 2 . correlation between manual and supervised learning approaches the kaplan - meier curves in figure 8 demonstrate the prognostic potential of our tacs-3 scoring system . tacs-3 negative patients showed significantly better disease - free and disease - specific survival compared to tacs-3 positive patients . in addition , cox proportional hazard regression showed significant correlation between our computationally generated tacs-3 scores and survival as listed in table 3 . although fiber features alone were correlated with survival , the highest correlation was observed when the tacs-3 scores were composed of multiple integrated fiber / epithelial features . this result shows that a multimodality imaging and analysis approach that combines features of not only collagen fibers , but both collagen fibers and cellular structures is most likely to succeed in predicting survival . km curves for disease free survival and disease specific survival showing the prognostic classification produced by our supervised learning based tumor associated collagen signatures-3 scoring approach univariate cox proportional hazard analysis results for various feature combinations table 4 lists the 14 most informative features in the tacs-3 scoring process ranked according to their weight produced by the linear svm algorithm . the svm weight was used to assess , which features were more or less informative in the classification . these features indicate the proximity between fibers and epithelial cell regions and were highly important in the tacs-3 classification . in addition , the difference in mean feature scores df = fp  fn for the training set is shown in table 4 for each of the ranked features . if df is > 0 , then the tacs-3 positive images had larger values for those features and if df is < 0 , then the tacs-3 negative images had larger values . for example , the density features resulted in lower df values in the tacs-3 positive cases indicating that the tacs-3 positive images had lower density collagen fibers . on the other hand , df was positive for the alignment features indicating that tacs-3 positive images tended to have more aligned fibers . interestingly , relative boundary angle was not as highly informative as many other features ; however , still was ranked within the top 14 of 27 features . feature ranking based on svm feature weight for a 16 patient ( 8 tacs-3 positive and 8 tacs-3 negative ) training set including the average feature value difference between the positive and negative training cases fp - fn the search for new prognostic and predictive breast cancer biomarkers is motivated by the need to improve patient outcome . a significant number of patients present with none of the currently available markers . in addition , survival and treatment response is often heterogeneous among patients within current biomarker classifications . the discovery and validation of new biomarkers will help to further improve breast cancer diagnosis and treatment planning . these new biomarkers need to be quantifiable , scalable and ideally correlate with both disease outcome and treatment specific response . the candidate biomarker we are focused on in this study ( tacs-3 ) measures collagen alignment relative to tumor - stromal boundaries and has been associated with progression in mouse models and has been shown to predict disease recurrence and survival in human patients . here , we demonstrate a protocol for using large field of view imaging techniques , image analysis and supervised learning to automate and quantify all of the steps in the process of tacs-3 scoring . these advances provide the tools for increasing the scale of tacs-3 investigations and applying tacs-3 scoring to cancers in other tissues such as ovarian and pancreatic cancer where these techniques could also be used to characterize other tacs both current and yet to be identified to see if they have research value in animal models or prognostic value in clinical specimens . tumor associated collagen signatures ( tacs ) analysis requires the simultaneous analysis of information about epithelial cells and extracellular collagen . the interactions between collagen and cells can only be assessed computationally if the cellular information is carefully registered with images of the collagen . we have therefore optimized our imaging system for highly automated capture of large fields of view , registered shg and bright field images of stained microscope slides with the purpose of analyzing collagen angle with respect to cell cluster boundaries . for this paper , we originally planned to use the same shg and bright field images captured by conklin et al . since these were already manually annotated . unfortunately , these images contained artifacts , which , although trivial for the human visual perception system to overcome , were extremely difficult for our computational systems to handle effectively . for example , shg images were originally captured in the backwards direction with elliptically polarized light , causing two artifacts . the first was simply a low signal to noise level due to few shg photons traveling in the backward direction from the thin tissue sections . the second artifact was observed as a larger relative shg signal from fibers in the direction parallel to the long axis of the laser polarization ellipse . artifacts in the bright field images included significant vignetting at field edges and low signal to noise due to short exposure times . these artifacts were easily hurdled by the human observers making tacs-3 assessments in a previous study . however , they are particularly difficult to handle by a computer vision based approach . we therefore decided to develop an optimized imaging system and protocol that would fix many of these artifacts and allow for more consistent automated imaging . similar image quality and consistency can be achieved with other shg microscopes including commercial systems with the appropriate hardware , but our analysis protocol did identify a necessary rigorous acquisition protocol that is best achieved with our new automated shg microscope . in general , the system should allow for fshg and bright field imaging with a field of view as large and flat as possible , numerical aperture of at least 0.75 , automated xyz motion control with appropriate position logging , circular polarization at the sample for shg imaging , autofocus and automated switching between shg and bright field imaging . our system of imaging and analysis to produce prognostic tacs-3 scores uses standard histopathology h and e slide preparations . the technique is therefore completely compatible with routine clinical protocols and is intended to augment currently available diagnostic tests . the current process requires no changes to current clinical protocol and the sample is returned to the clinician unmodified . we present a system that uses shg imaging to capture collagen fiber images ; however , wide field polarization sensitive techniques such as lc - polscope or picrosirius red staining might be used to alternatively capture images of collagen fibers . one advantage of using shg is that it does not require additional stains and can capture three - dimensional fiber information in thick , unstained tissue samples . unfortunately , when imaging in thick unstained tissue , the identification of epithelial regions can be difficult ; however , techniques using autofluoresence and fluorescent lifetime imaging have been shown to be capable of this task . as implemented here , our tacs-3 scoring algorithm is necessarily two - dimensional , since we are relying on h&e stained slides for our epithelial cell information . however , fiber extraction , epithelial region segmentation and relative angle measurements can be extended to three - dimensions without significant alteration of our general protocol . in addition , although our current tacs-3 scoring protocol is able to process standard h&e stains , staining for epithelial cells , with , for example , pan - cytokeratin conjugated stains , may simplify and improve epithelial cell segmentation . future methods may also be adapted to segment clusters of fibroblasts , macrophages and other stromal cells , whose proximity and relative morphological structure with respect to surrounding collagen fibers may further improve correlation with survival or metastatic potential . collagen alignment related image features are interesting not only because they have been shown to be prognostic , but because they have been shown to be directly linked to cancer biology . researchers have found that cells are more likely to invade along parallel , aligned collagen fibers , features that are directly being measured by our system . access to the breadth of fiber data available with our techniques could lead to advances in our understanding of these biological phenomena . relevant feature sets are not always available with other machine vision systems developed for biological image classification . for example , although wndchrm is an extremely powerful image classification tool , informative image features often do not relate to the biology at hand . in the case of our tacs-3 analysis system , biological observations have driven the image analysis model ; therefore , features are more easily linked back to biological functions potentially revealing new insights . high mammographic density is one of the largest risk factors for the development of breast cancer and has been associated with increased epithelial cellularity and increased collagen density . increased collagen density has been observed to promote tumor progression in a mouse tumor model and in node positive breast cancer leading one to potentially conclude that collagen density causes elevated risk . however , maller et al . observed that high density , nonfibrillar collagen protected against tumor progression and alternatively , that linearized collagen fibers induced invasive cellular behavior . in agreement with these recent findings , we observe here that tacs-3 fibers are more commonly present in regions of lower fiber density and are more likely to be thinner , more linearized fibers . thick , curvy , and denser collagen fibers are unlikely to contain tacs-3 fibers and are observed to be associated with a better prognosis . these observations support the hypothesis that collagen fiber shape and organization is a key aspect of the invasive extracellular matrix ( ecm ) phenotype . the imaging instrumentation presented here consists of a relatively compact and highly automated multiphoton microscope with an integrated bright field slide scanner . the system has been optimized to capture registered whole slide images of both bright field and shg images of histopathology slides by imaging each small  20 field of view and automatically aligning and stitching each image together . capturing large fields of view in this manner allows for a more thorough and consistent data collection potentially reducing sampling bias and supporting pipelined computational image analysis . the registration of cellular with extracellular collagen information provided by our system allows for the quantitative analysis of key relative structural features between collagen fibers and cancer cell clusters . in addition to shg , our multiphoton system is capable of imaging other endogenous fluorophores such as nicotinamide adenine dinucleotide ( nadh ) or flavin adenine dinucleotide ( fad ) as well as any of the routine exogenous multiphoton probes used to stain tissue . we present an imaging and analysis protocol that uses high content imaging techniques coupled with supervised learning to perform semi - automated tacs-3 scoring of slide mounted biopsy samples . we apply our technique to a previously annotated tma containing tissue from 207 patients with invasive breast cancer . the resulting scores are shown to positively correlate with manual annotations and to predict patient outcome with good statistical significance . future work will attempt to validate this technique on larger cohorts of breast cancer patients , to study ecm targeted drug responses in animal models , and to study collagen alignment in other cancers . as well , future work will focus on improving the clinical application of these techniques so they can be run by untrained clinical personnel and be run at the time of acquisition to find rois and tacs within those regions automatically . together with more automation , tacs screening has great potential as a clinical diagnostic tool that can provide relevant prognostic information from large numbers of tissue samples .
background : mortality in cancer patients is directly attributable to the ability of cancer cells to metastasize to distant sites from the primary tumor . 
 this migration of tumor cells begins with a remodeling of the local tumor microenvironment , including changes to the extracellular matrix and the recruitment of stromal cells , both of which facilitate invasion of tumor cells into the bloodstream . in breast cancer , it has been proposed that the alignment of collagen fibers surrounding tumor epithelial cells can serve as a quantitative image - based biomarker for survival of invasive ductal carcinoma patients . 
 specific types of collagen alignment have been identified for their prognostic value and now these tumor associated collagen signatures ( tacs ) are central to several clinical specimen imaging trials . here , we implement the semi - automated acquisition and analysis of this tacs candidate biomarker and demonstrate a protocol that will allow consistent scoring to be performed throughout large patient cohorts.methods:using large field of view high resolution microscopy techniques , image processing and supervised learning methods , we are able to quantify and score features of collagen fiber alignment with respect to adjacent tumor - stromal boundaries.results:our semi - automated technique produced scores that have statistically significant correlation with scores generated by a panel of three human observers . 
 in addition , our system generated classification scores that accurately predicted survival in a cohort of 196 breast cancer patients . 
 feature rank analysis reveals that tacs positive fibers are more well - aligned with each other , are of generally lower density , and terminate within or near groups of epithelial cells at larger angles of interaction.conclusion:these results demonstrate the utility of a supervised learning protocol for streamlining the analysis of collagen alignment with respect to tumor stromal boundaries .