Quantitative analysis of second harmonic generated images of collagen fibers: a review

The human body is a complex structure. Its strength is ensured by the collagen protein which exists under the form of fibers. The quantitative analysis of these fibers in biological tissues can be very interesting to establish a relationship between the microstructure and their functions. This analysis is usually performed using two-photon microscopy and second harmonic generated (SHG) images. Lately, more and more researchers focused on the use of SHG images since it is a non-invasive technique and allows the capture of collagen fibers only. Many image-processing techniques can be used to extract quantitative information from those images such as fiber orientations, dimensions, and density. Therefore, accurate measure extraction depends mainly on the used image processing methods and, thus, it is necessary to know what processing technique to use. The main purpose of this article is to exhibit the most used techniques in collagen fiber quantitative analysis then categorize them according to the information to extract. A comparison of three most used methods in fiber orientation’s estimation is carried out. Despite the considerable number of papers aiming to quantitatively analyze collagen fibers from SHG images, two main aspects were not deeply covered. First, the use of deep learning algorithms is still limited even for segmentation and denoizing applications. Second, most of the studies processed in this review focused on two-dimensional SHG images and did not take into consideration collagen fibers as a three-dimensional volume.


Introduction
Collagen is the most abundant protein in the human body and in mammals, in general. This protein is what holds the body together since it ensures the strength and elasticity of the body's connective tissues. It can be divided into different types. Around 80 to 90% of the collagen in the human body are collagen types I (skin, tendons, bones, ligaments), II (cartilages), and III (scar tissue, muscles, vessel walls) (Lodish et al. 2000). It is mostly found in the form of fibers, but it is important to recall that the structures which can be seen in microscopic images depends on the imaging scale and what is observed: one can see bundles composed of collagen fibrils at a 20-µm scale or single collagen fibrils at a 100-nm scale. The fundamental structural unit of these fibrils is a triple helix with a length of 300 nm and a diameter of 1.5 nm (Lodish et al. 2000). For simplicity, the bundles are also often called collagen fibers, and their typical dimensions are in the order of 1 to a few tens of micrometers in diameter and several hundreds of micrometers in length. The study of these fibers, which are essential to the proper function of tissues, is fundamental in understanding the etiology of pathologies, their evolution, and in improving their clinical diagnosis and management. It is a multidisciplinary field involving mechanics, image processing, chemistry, biology, etc. Relevant studies targeting this aspect require both suitable imaging techniques and reliable image analysis methods.
Early studies conducted on this protein led to the characterization at the microscopic scale (around 30 µm) of many tissues composed essentially of collagen by the mean of histology. This technique consists in studying the microscopic structure of biological tissues and the relations between individual elements (Lowe et al. 2015). It involves a chemically destructive process, which can only be performed on ex vivo samples, and eventually a slicing step, which prevents from 3D observations. The process can also have an impact on the microstructure of the considered specimen. For example, freezing of the tissue is used for visualization purposes, which may destroy some of its components. Its use remains, however, a standard for pathology diagnoses in clinics. Hence, the study of the microstructure evolution of a biological system ex vivo under a mechanical load is impossible with histology. In order to deal with this issue, new imaging techniques have been tested and proved their efficiency. Among those imaging modalities, one can cite first scanning electron microscopy (SEM) (Prado et al. 2003;Orberg et al. 1982). It allows obtaining images with a resolution of 1-20 nm, but it lacks in-depth signal: only the peripheral surface can be imaged accurately. As a remedy, other techniques would ideally ensure a precise quantification of the collagen fibers within the volume of the material as oriented structures in space in an adapted scale (1-100 µm). For instance, X-ray computed tomography (XRCT) and X-ray micro-tomography are well suited for quantifying collagen fiber architecture because they allow capturing their microstructure through larger fields of view (up to 1.7 mm × 1.7 mm) as compared to other microscopy techniques (Bailly et al. 2018;Disney et al. 2017;Walton et al. 2015). In addition, it offers a resolution of 20-100 µm, though a compromise between resolution and field of view must be made. However, the addition of X-ray contrast agents may change the behavior of the specimen components, restraining their use. Optical coherence tomography (OCT) (Fujimoto et al. 2000) has been used as an alternative (Babalola et al. 2014;Ugryumova et al. 2009). Just like XRCT, OCT provides a resolution of 1-15 µm, but it does not allow capturing individual components of a specimen. This makes quantitative analysis hard to achieve such as for aortic ostial lesions where it is not possible to clear the blood at the entrance to neighboring arteries. It is dependent on the considered biological tissue scattering and absorption. Yet, it is possible to use optical clearing agents to reduce light scattering but it can have an impact on the tissue structure. Recently, a strong interest was shown toward fluorescence microscopy, which requires the use of stains, but has limited physico-chemical modification of biological tissues. For instance, confocal microscopy was often used (Wu et al. 2003;Stein et al. 2008) because of its resolution of around 160 nm and its capability to capture images through the specimen depth. Later, the emergence of powerful lasers enabled multi-photon microscopy. This imaging technique, with (Chen et al. 2012;Polzer et al. 2013;Yeh et al. 2002) or without polarizer (Ayyalasomayajula et al. 2019;Cavinato et al. 2017) does not harm the sample because it is less exposed to the laser. It offers a scale for representation of the order of a micrometer and a resolution up to 150-200 nm. Besides, it allows imaging deeper into the sample and thus collecting more images in the depth (up to 500 µm (Yamada et al. 2014)). Additionally, collagen fibers react to multi-photon laser by generating second harmonics across the spectral region between 400 and 500 nm (Theodossiou et al. 2006). This property, called second harmonic generation (SHG), is an asset to capture images of collagen fiber only, as this signal is specific and can be separated from other signals.
The cited imaging modalities introduced some improvement on how to capture sufficiently good images to extract information related to the structure and the function of collagen fibers. Studying the organization of the collagen fibers is of interest in biomedical research since it allows diagnosing fibrosis (Campagnola 2011;Strupler et al. 2007) or analyzing their interaction with cancerous cells (Bredfeldt et al. 2014). Moreover, in Brown et al. (2003), the authors tried to quantify the dynamics of collagen modification in tumors in vivo after pharmacologic intervention. Besides, other researchers focused on the quantification of fiber orientation in order to study structure-to-function relationships such as in pressurized vessels (Ayyalasomayajula et al. 2019;Cavinato et al. 2017;Schriefl et al. 2013) or waviness and density in order to identify the impact of sample aging (Sugita and Matsumoto 2017; Wu et al. 2016).
However, this type of quantification is complicated and requires dedicated methods. To date, the conducted research in this field succeeded in characterizing collagen structures only in 2D planes although those tissues are three dimensional.
In the present review, we focus on the analysis of those SHG images of collagen fibers and the different techniques developed to extract information aimed at characterizing those fiber networks. We emphasize on the accurate quantification of several quantities such as the fiber orientation, waviness, and dimensions in addition to the collagen density. To this aim, the paper is organized as follow: the first section will explain the physics behind SHG imaging. The second section will focus on the different image processing techniques used to quantify collagen fibers in SHG image. The third section will exhibit what has been done in the literature to extract the collagen fiber geometric, composition, and morphological information from SHG images. This review will end with a comparison between of three common methods used to extract fiber orientation from SHG images of collagen networks.

Study design
Articles that could answer questions relating to the quantitative analysis collagen fibers were carefully selected. Only papers dealing with biological tissues and second harmonic generation acquisition technique were taken into consideration. Papers dealing with polarized SHG and papers with only abstracts were excluded. The search was performed using some of the most popular digital repositories (Google Scholar, PubMed, Springer, Science Direct and Optica publishing group) by using keywords "SHG," "collagen," and "quantitative analysis." The collected articles were divided into three categories with respect to the quantitative information extracted (the geometry and morphology of collagen fibers and the composition of the considered specimen in collagen). Then, these articles were categorized following the techniques used to improve SHG images or to extract relevant information. The most frequent techniques were selected to be cited and explained in this review.

SHG acquisition technique
Optical harmonics were first discovered in the 1960s when the high-intensity pulsed lasers have been invented. Franken et al. (1961) observed second harmonic generation (SHG) in crystalline quartz by using a Q-switched ruby laser. It became a very used method to characterize the second-order nonlinear optical (NLO) response of emerging materials, especially organic NLO materials. It led to an increase of its use in different fields such as biomedical research (Singer and Wu 2013). SHG imaging in biology was reported by Freund et al. (1986) when he tried to characterize the polarity of collagen fibers in a rat-tail tendon. Campagnola et al. (2002) reported a more recent practical implementation in where the authors succeeded in imaging biological tissue at high resolution and fast acquisition rate. In order to collect SHG images of collagen in biological tissues, a two-photon light microscopy has been developed. This imaging technique is based on the excitation of a molecule to a virtual state by two photons which are then converted into a single photon of same total energy at double frequency, without absorption or re-emission of photons. For two-photon excitation, photons in the infrared spectral range are used under highly intense laser illumination (for example, Ti:sapphire lasers). Infrared photons are chosen because of their low energy. When the energy gap between the ground state and the excited state is smaller than the sum of the energy of the two photons, the nonlinear process can occur. In this case, the probability that a fluorescent molecule absorbs simultaneously two infrared photons is a quadratic function of the excitation radiance (So 2002).
The possibility to take microscopic images in three dimensions (i.e., depth discrimination) is considered one of the most interesting properties of two-photon microscopes. It originates from the almost absence of out-of-focus light resulted by the reflection. Eighty percent or more of the total florescence signal may be cramped in a region of 1 µm thickness around the focal plane of a two-photon excitation (So 2002). Notable two-photon excitation occurs where the photon density is high. It corresponds to the focal volume of the microscope, which can be as small as ~ 0.1 µm 3 (Zipfel et al. 2003). The laser needs to scan the entire specimen in the three dimensions to generate a 3D image, which may involve relatively long acquisition times according to the volume size and acquisition parameters. Finally, the use of two lowenergy photons limits the risk of photo-damage of the sample (Svoboda and Yasuda 2006). In addition, it maximizes the probability of detecting photons per excitation event in the right spot and, thus, minimizes photo-bleaching (when the molecule loses its ability to fluoresce) and photo-toxicity.

Method description
We here provide a succinct description of the methods allowing extracting quantitative information from SHG images of collagen fibers. We first introduce some common pre-processing techniques. Then, we focus on the different image transformations that may be used to analyze the SHG images. Finally, we highlight the methods useful to extract and select information from SHG images. Fig. 1 displays the methods that will be covered in this section.

Pre-processing
In order to process an SHG image or stack and to extract as much accurate information as possible, it is important to remove the noise. SHG images of tissues present usually Poisson noise because of their poor signal-to-noise ratio (Bredfeldt et al. 2014). Common methods in pre-processing SHG images are introduced in the following and are summarized up in Table 1.

Median filter
It is simple and widely used filter to reduce noise in images and to smooth them. Median filter (Huang et al. 1979) is a nonlinear smoothing filter. The value of each pixel in the image is replaced by the median value of pixel's intensity in a previously defined neighborhood of size m * m.
Median filters work well for removing random salt and pepper noises (Gonzalez and Woods 2018). However, this kind of filters does not allow the suppression of Gaussian noise (Ohki et al. 1995) which can be dealt with through deconvolution. They do not reduce the difference in brightness of images and, hence, preserve edges. However, when the signal-to-noise ratio of the image is small, or the neighborhood is too large, median filters tend to delete useful information and produce false noise edges.

Contrast enhancement
The recognition of image features depends on the image contrast. However, the contrast can be distorted by the imaging system because of poor illumination conditions. For this purpose, histogram equalization is widely used. A well-acquired gray-scale image should cover black and white pixels. It is also better that the image's shades are evenly distributed (i.e., the image histogram is uniform). Many contrast transforms can be used for this purpose such as histogram equalization, adaptive histogram equalization, and contrast limited adaptive histogram equalization (CLAHE) (Mustafa and Abdul Kader 2018). Here, we will focus on the CLAHE algorithm, which is very used on SHG images (Koch et al. 2014;Hu et al. 2012).
CLAHE is a variant of adaptive histogram equalization (Pizer et al. 1987). It consists in computing histograms of distinct regions and using them to redistribute the pixel intensity values of the image. The difference between CLAHE and other adaptive histogram equalization algorithms is that it clips the histogram at a pre-defined value (i.e., if a histogram bin is higher than the clip limit, those pixels are clipped and uniformly shared with other bins before proceeding to the histogram equalization). It operates on small regions of the image called tiles. To remove the artificial boundaries between the different tiles, bi-linear interpolation is used.
CLAHE is a good technique to improve local contrast and to enhance edges. Compared to other adaptive histogram equalization techniques, it limits the noise amplification. However, if the image is too noisy, a phenomenon of noise amplification may occur. The combination of CLAHE, median filter, and edge sharpener (such as high-pass filters) can be successful to maintain the image high spatial frequency content.

Directional filters
When there is a need to study oriented features in an image, directional filters can be used (Bamberger and Smith 1992). They consist in a filter bank containing lines in different directions. They can be used to detect edges or to identify objects orientations. Those filters have wedge-shaped passband spectral regions, and are therefore usually referred to as wedge or fan filters (Simoncelli and Farid 1996). When the orientation of the wedge is known, the determination of objects direction is straightforward.
Wedge filters are an easy to implement and efficient tool to study oriented objects in images. However, in order to have a fine description of the image features, it is necessary to have thin wedges in as many directions as possible which will add more computational costs.

Gradient
The gradient vector is a fundamental approach for finding extrema of a continuous and smooth function in space (Hyvärinen et al. 2009). The gradient is defined as the partial derivatives of a function with respect to all its components. For 2D images, the gradient is usually achieved by the convolution of the image by a couple of filters based on the Sobel filter (Sobel 1968) or the Prewitt operator (Prewitt 1970).
The fact that small displacements are taken into consideration to compute the gradient allows capturing as much details as possible in any direction. The gradient works fine with clear images without much noise. However, if the considered image is noisy, the gradient will not bring any useful information.

Frangi filter
The Frangi filter (Frangi et al. 1998) was first developed to be a vessel enhancement filter. However, it was used to detect both vessel-like and tube-like structures in images. Because of the collagen fiber morphology, which can be assimilated to tubes, the Frangi filter was used to extract the fibers from SHG images.
The Frangi filter is based on the computation of the image's Hessian matrix. In the proposed framework, the derivative of an image corresponds to its convolution with derivatives of Gaussians. The second derivative of a Gaussian kernel shape allows to measure the contrast between the region in and out of a range (− s, s), s being the standard deviation of the Gaussian. Through the analysis of the eigenvalues of the image's Hessian matrix, it is possible to extract the direction of the smallest vessel's curvature which corresponds to the main directions in which the local second-order structure of the image can be decomposed (Frangi et al. 1998). The eigenvalue decomposition gives three orthonormal directions, which allow describing vessels in images.
The use of multiscale Frangi filter, through the analysis of the eigenvalues of the Hessian matrix, makes it possible to capture the smallest details of an image and, thus, avoid the application of different filters of different sizes. However, the Frangi filter may not take into consideration any object in the image, which does not have a circular cross-section.

Image transformations
A strong interest was shown to signal decomposition because of the uneven distribution of signal energy in the frequency domain. It consists in dividing the signal spectrum into its sub-spectra, which are then treated individually (Akansu 2001). Signal decomposition was used for many applications such as compression and feature extraction. For image analysis, and particularly for studying the collagen fibers in SHG images, several image decomposition methods have been used: the fast Fourier transform (FFT), the wavelet transform (WT), the radon transform (RT), and the Hough transform (HT). Table 2 sums up these methods.

Fast Fourier transform (FFT)
The FFT is an efficient method to study the spatial frequency distribution of the pixels in an image. The Fourier transform (FT) was initially used to characterize linear systems and to identify their frequency components that make a continuous waveform (Bergland 1969). Images are processed using the discrete Fourier transform (DFT). The DFT coefficients can be computed by the FFT. This transform is a computationally cheap and fast algorithm originally introduced by Cooley et al. (1969). Different approaches can be chosen to compute the FFT (Rader 1968;Bluestein 1970;Bruun 1978;Rivard 1977). The 2D DFT is a two-time process (Rivard 1977). It consists in combining vertical and horizontal 1D DFT of an array into one 2D transform that makes sense. First, a 1D DFT over horizontal lines of an image is performed. Then, a 1D DFT over vertical lines is applied on the result of the previous operation.
The 2D FFT is an efficient operator to characterize an image and to capture the variation of its texture. However, the space notion is lost when the transition from space to frequency is done. In fact, the 2D Fourier transform gives information of the global contents and changes in frequency without knowledge on the section of the image that corresponds to it. Besides, Fourier transform may not work accurately to reconstruct an image, which is highly non-smooth (Jaffar Iqbal Barbhuiya and Hemachandran 2013).

Wavelet transform (WT)
Wavelet methods have become a widely used tool in image processing during the last 20 years. This is due to their ability to analyze non-stationary structures and characterize local properties. An image is mapped to a phase space, which is  (Frangi et al. 1998) -It only takes into consideration objects with circular cross-section (Frangi et al. 1998) parameterized by a scale/size/resolution and a time/space parameters. Wavelet transform is an alternative to the Fourier transform which characterizes the image in a time/space frequency space (Dahlke et al. 2008). It was first introduced by Grossmann and Morlet (1984) as an elegant multi-resolution signal processing tool thanks to its ability to naturally vary the time-frequency resolution (Akansu 2001). It is a mathematical function of zero average used to divide a function into components at different scales. Each scale is computed using a specific wavelet generated from an initial function named mother wavelet by dilation and translation. The dilation allows carrying out a multi-scale analysis and enables to capture small details. It is possible to perform a wavelet decomposition of an image (in 2D or even for higher dimensions) in order to compress the data or to obtain a vector of features that characterizes the data in a basis of wavelet. It is helpful to capture the orientation changes in an image. For this matter, we need to perform a 2D discrete wavelet transform (DWT). As for the 2D FFT, it can be generated using the horizontal and vertical 1D DWT.
The main advantage of the wavelet transform is that it provides a localization in both space and frequency domains. The wavelet transform allows capturing small and coarse details. Indeed, wavelet transforms over-perform traditional Fourier transforms in representing functions with sharp peak discontinuities and in correctly decomposing and reconstructing non-stationary, non-periodic, and finite signals (Jaffar Iqbal Barbhuiya and Hemachandran 2013). It can also be used to detect discontinuities and irregularities in signals. However, this technique is computationally expensive for fine decomposition. The choice of the mother wavelet and the number of decompositions can highly influence the result.

Radon transform (RT)
The radon transform is a mathematical transformation based on projections, which is the basis of computed tomography (CT). We can also use it to detect edges. The radon transform consists in performing different projections of an image according to different angles. The resulting projection corresponds to the integral of the line integral (i.e., the sum of the pixel intensities in every direction) (Deans 2007). In other terms, the RT maps an image from Cartesian coordinates to polar ones. This transform can also be applied to 3D images. In this case, the integral is taken over planes. The RT data is usually referred to as sinograms.
The use of a FFT gives qualitative information about fiber orientations. To deal with this issue, it is possible to apply an RT on the result of the FFT. Since it is based on projections, it gives quantitative information for each considered angle. It is important to have a sufficient number of angles to get accurate results in detecting and extracting the fiber orientations.

Hough transform (HT)
The Hough transform (Hough 1962) was first introduced to detect lines in images. This algorithm was then simplified by Duda and Hart (1972) and generalized to detect circles and curves. The original HT algorithm assumes that every line in an image can be represented by a unique couple (slope, intercept). Duda and Hart (1972) changed this representation by the couple (angle, distance), where the angle and the distance correspond to the polar coordinates of a considered line in the image (the distance being the distance between the image origin and its projection on the line). A matrix called accumulator is created where its axes correspond to the parameters characterizing the line. Thus, for each pixel of the image, the accumulator is incremented for all possible lines passing through that pixel. The presence of an edge -The result highly depends on the choice of the mother wavelet -It is computationally expensive for fine decomposition RT Projection data -It gives information with respect to the angle of projection (Deans 2007) -It depends on the chosen number of angles (Deans 2007) -It is computationally expensive for fine analysis HT Polar map of the image -It corrects properly the detected edges (Leavers 1992) -It can be used to estimate object orientations -It works better on detected edges -It does not distinguish between objects if they are aligned corresponds to a high value position in the accumulator (Leavers 1992). A reconstruction of the initial image is possible by retrieving the parameters corresponding to the peaks in the accumulator. The HT gives good results when applied on an image where the edges were already detected. It works fine with noisy data. This method allows reconstructing an edge if it is discontinuous when performing the edge detection algorithm. Therefore, the application of the HT needs a prior step to detect the edges. For linear objects, the HT is a good method to detect edge orientation directly from the accumulator matrix. However, its effectiveness depends on the considered image: if two objects are aligned in an image, the HT will exhibit them as one.

Information selection and extraction
After the pre-processing of an image, the SHG image analysis needs to extract as much valuable information as possible. For this matter, it is possible to extract this information through a spatial characterization or a statistical one. For both types of characterizations, many methods can be used. Some of them are detailed hereafter and are summed up in Table 3.

Spatial information selection
To analyze an image, it is important to consider the spatial distribution of the pixel intensity. This is possible through several techniques: (i) segmentation, which transforms an SHG image into a binary image where only the collagen fibers are represented; (ii) skeletonization, which determines the center line of the collagen fibers in the SHG images and thus, allows extracting geometrical information about the fibers.
Pixel-based segmentation This type of segmentation aims to gather pixels corresponding to an object and mark them. It is based on their intensity similarity and spatial proximity. The (automatic) thresholding segmentation is the easiest method for image segmentation. Otsu thresholding algorithm (1979) is the most used one, especially on SHG images, because of its simplicity in addition to the fact that it works particularly well when the considered image contains two classes (an object and the background). Its principle is to find the threshold that maximizes the interclass variance of a two-classes histogram. In addition to this method, several other approaches exist to compute the threshold such as entropy-based thresholding (Khattak et al. 2015;Luthon et al. 2004), minimum error thresholding (Kittler and Illingworth 1986), moment-preserving thresholding (Tsai 1985), fuzzy set thresholding (Tizhoosh 2005), etc.
Thresholding decomposes the image gray scale information with respect to gray level of targeted objects. There are two types of thresholding segmentation: global and local. The global threshold looks at the global picture: it divides the image into two regions (background and target) regardless of the positions of objects. On the contrary, local thresholding looks for a threshold in a neighborhood around any pixel of the image.
The main advantages of thresholding techniques are their simplicity and their fast computation. This type of segmentation works well when the image's histogram presents two or more peaks. However, it is highly sensitive to the tackled problem and is specific to the considered image. In addition, it takes only into consideration the intensity of the pixel/voxel and not its spatial information, which makes this method highly sensitive to noise (Yuheng and Hao 2017). In fact, small areas or isolated pixels can be classified as independent regions even though they represent noise or belong to another region. Besides, to segment SHG stacks where the pixels intensity decreases with depth, it is complicated to find a threshold that takes into consideration that phenomenon.

Region-based segmentation
Unlike pixel-based segmentation which classifies a pixel-based on its intensity value without taking into consideration the spatial context, region-based segmentation looks for pixels having similar features. Several techniques belong to this category such as region growing algorithm (Adams and Bischof 1994;Mancas et al. 2006), split and merge algorithm (Damiand and Resch 2003;Chaudhuri and Agrawal 2010), and clustering (Thilagamani and Shanthi 2011). Our interest is paid to the region-growing algorithm since it has been used in quantifying SHG images of collagen. First, the user selects initial seed points to be in a region. Then the algorithm checks iteratively if the adjacent pixels should be added to the region according to one or several of available criteria (gray scale texture, intensity, color, etc.) (Yuheng and Hao 2017).
Region-based segmentation allows partitioning the image into sub-regions. However, those methods depend on the choice of seed points and do not work properly on nonsmoothly varying regions. Besides, a threshold is needed as a criterion to construct the regions; thus, its choice is important. Finally, it is a local technique with no global view on the image and it is sensitive to noise, which may lead to an over-segmentation.
Edge-based segmentation An important feature carrying information about objects is their borders, i.e., the discontinuities in the pixels' intensity. To detect the gray level discontinuities, the most common approach is based on detecting edges, which represents a set of connected pixels forming a boundary between two regions (Gonzalez and Woods 2018). There is a gap between the pixel values of two adjacent regions. Those discontinuities can be either step edges or line edges.
Step edges are characterized by the sudden change in the pixel intensity from a region to another. Line edges correspond to a sudden change of the pixel values followed by another sudden change to return to the initial value within a short distance (Senthilkumaran and Rajesh 2009). However, in real images, it is impossible to find those types of edges because of the smoothing introduced by the optical systems or by the low-frequency components of images. One can find ramp edges instead of step edges and roof edges instead of line edges, where the pixel intensity change occurs over a finite distance. Such gaps can be detected with the help of differential operators such as the Sobel operator (Sobel 1968), the Laplacian and the Laplacian of Gaussians (also called Marr-Hildreth operator) (Acharya and Ray 2005), the Prewitt operator (Prewitt 1970), or the Kirsch operator (Kirsch 1971). More sophisticated techniques such as the Hough transform were also used to determine image edges (Hough 1962). Once the edges are detected, mathematical morphology operators (erosion, dilation, opening, closing, etc.) are used to fill the targeted regions and, thus, segment the image.
Edge-based segmentation is a high-level segmentation approach similar to the way humans perceive an image. It works well on images with high contrast. However, it is highly sensitive to noise. It is centered on local information and does not take into consideration the global view. In addition, it does not work well to detect corners and when the contrast is low.

Fast marching method (FMM)
The fast-marching algorithm (Malladi and Sethian 1996) allows to track object boundary. It was initially developed to follow an interface or contour propagating under a speed function F and was then used in medical applications (Cardinal 2010). The FMM is a discretized and computationally optimized version of the level set method (Osher and Sethian 1988). It aims at spreading an initial surface until it covers the entire surface of interest (the collagen fibers in our case) by solving the Eikonal equation. It is based on computing a distance map between the initial surface and its surroundings. The surrounding points are divided into three regions: the accepted points, the narrow band, and the far region. Initially, the accepted point region is the initial surface. The narrow band constitutes the closest pixel to the initial front. The far region is what is left of the image. The Eikonal equation is solved on the edge points of the initial surface. The points that satisfy this equation are then added to the initial surface and the same steps are applied again until there are no more points that may be added to the accepted point set.
The algorithm gives good results when the image is very distinct from its background. Besides, the use of such algorithm does not need a prior setting of the parametric representation of the surface contour to be followed: this technique is robust with respect of the topology to be analyzed. However, it relies entirely on a physical interpretation of the problem characterized by the isotropic front propagation of the initial surface (Cristiani 2009). Besides, the use of the first-order neighbors (only four neighbors) introduces errors in the computation of the travel time from a point to another.

CT-FIRE
The CT-FIRE is an algorithm introduced by Bredfeldt et al. (2014) that enables the extraction of fibers through their skeletons. It was developed to extract collagen fibers from SHG images in order to estimate their orientation and geometric information.
This algorithm is based on two steps. The first one is a filtering using curvelet transform (CT). Curvelet filters were introduced by Starck et al. (2002) in order to overcome the limitation of highlighting lines and edges. The curvelet transform is a wavelet transform except that instead of the wavelets we use curved functions called curvelets. The second step consists in applying the fiber extraction algorithm FIRE developed by Stein et al. (2008) It describes the fibers as a set of n vertices and p paths. Every path corresponds to a fiber characterized by k vertex identifiers p i = (n i 1 , n i 2, …, n i k ). The image is smoothed using a Gaussian filter before segmenting it through thresholding. Then, for each pixel of the segmented image, the Euclidean distance map is computed. This map is used to identify the centerlines of the fibers. Once the centerlines identified, short non-relevant fibers are deleted and close fibers are connected.
The CT step introduced in the CT-FIRE algorithm improved the result of the fiber extraction compared to the classic FIRE algorithm. It provides better results when the collagen fibers are densely packed. However, for highly noisy images, other pre-processing techniques may be needed before applying the CT-FIRE algorithm. It also does not work well on images where the fibers are wavy and present many intersections.

Statistical feature extraction
The analysis of an image texture covers the region-specific identification of higher-order properties which are hard to detect visually. Texture analysis leads to the definition of statistically uniform regions of an image based on the intensity distribution (Dudenkova et al. 2019). Statistical approaches that have been used to analyze SHG collagen images can be divided into three categories: firstorder statistics, second-order statistics, and directional statistics.

First-order statistics (FOS)
First-order statistics estimates parameters derived directly from the image statistics. They are often used to simply describe the image intensity distribution. However, they ignore the spatial correlations between the pixels of the image. In other terms, FOS describes the probability to observe a pixel having a certain intensity in any position in the image. In more details: • The intensity distribution histogram is a representation of the number of pixels in an image with respect to their values. It is a useful tool to detect saturation effects in an image (i.e., presence of pixels with maximum intensity), to deduce the brightness (the image is bright if the histogram values are more concentrated around high values), and to check the contrast (if the values of the histogram are spread out without a noticeable peak).
• The mean calculated from the pixels' intensity or from the probability distribution of the pixels' intensity is used to evaluate the presence of one texture in the image. • The standard deviation captures how the pixels are spread out with respect to their intensity. • The skewness evaluates the histogram's lack of symmetry and allows characterizing the slope of the image histogram with respect to the central line. The skewness of a normal distribution is equal to zero. A negative (resp. positive) skewness denotes an image for which the majority of pixels have values smaller (resp. greater) than the mean value. • The kurtosis describes how much a distribution is concentrated around a peak (the mean) and allows evaluating the efficiency of a denoizing algorithm.
FOS are easy and fast to calculate. However, their interpretation is not always simple. They give global information and cannot be used to quantify local information (unless the initial image is divided into several ROIs).

Second-order statistics (SOS)
Second-order statistics estimate parameters from the matrix generated by performing a correlation between the image pixels. It studies, in particular, the topology of one region compared to the image. Here we talk about texture analysis. This technique is usually used to describe and characterize a local area in an image through the use of gray level co-occurrence matrix (GLCM) and some statistics (Haralick et al. 1973).
• The GLCM evaluates the spatial relationships between the values of the pixel intensity. It is a squared matrix of dimension equal to the number of gray levels in the considered image (for example, 256 for 8-bit images). The parameters that will be presented subsequently (Iqbal et al. 2021) can be calculated from the initial image but they are more relevant when they are performed on the GLCM. • The energy (also called uniformity) allows to evaluate the uniformity of the image. • The inverse difference moment (IDM) measures the local homogeneity of an image. When the IDM value increases, it means that the incidence of pixels' pairs cooccurrence is enhanced which means that IDM is high when the image is homogeneous. • The inertia (also called contrast) allows studying local variations in an image. It is highly sensitive to large differences in the GLCM values and has a strong correlation with the lowest and highest values in a ROI. • The correlation characterizes the gray-level linear dependency on specified pixels on an image (i.e., the repetitive nature of the texture element position).
• The entropy focuses on the randomness of regions in an image with respect to its neighborhood in terms of intensity distribution. Low entropy values correspond to a uniform and homogeneous image Texture analysis can be performed on an entire image, but it is more interesting on a localized area to capture morphological changes. This technique allows seeing morphological changes of the collagen structure (for example, to make a comparison between a benign and a malignant tumor), but it does not give information about their geometric and composition information.
Directional statistics Directional statistics focuses on observations that have directions. These observations usually lie whether on the circumference of a circle (circular statistics) or on the surface of a sphere or a hypersphere (spherical statistics) (Ley and Verdebout 2017). Statistical analysis of directional data became more used after Fisher's paper (1953) where he explained the need to consider the curved nature of the sample space. Several directional distributions emanated from Fisher's contribution. They are based on the extension of classical concepts from multivariate analysis (e.g., point estimation, regression, multi-sample testing procedure) to directional setting (Pewsey and Garcίa-Portugués 2020; Mardia et al. 2008;Mardia and Jupp 2000). In the following, we will focus on the Von Mises distribution which has been used to extract quantitative information from SHG images of collagen fibers.
The Von Mises distribution is considered a flexible circular distribution. It is useful for a circle from a statistical inference point of view (Mardia and Jupp 2000). It represents the maximum entropy distribution for circular data when the first circular moment real and imaginary parts are specified. It is characterized by two parameters, a location parameter µ ∊ [− π, π] and a concentration parameter κ. κ is positive and it allows to regulate the concentration of the distribution around µ. This distribution was later generalized to higher dimensions by Von Mises and Fisher and, thus, was named von Mises-Fisher distribution. The Von Mises distribution can also be referred to as the circular normal distribution. To characterize collagen in SHG images, it is possible to evaluate the fiber dispersion and its diameter by fitting a Von Mises distribution.
It is a good tool to study 3D images because it can be generalized to high dimensions without using many parameters. However, for SHG images of collagen, this method assumes that all the fibers belong to a single family (i.e., having the same orientation).

Quantities to measure and associated questions
In order to analyze and understand how collagen fibers behave when they are under a mechanical load, it is necessary to quantify them using some relevant information. For example, for arteries, the quantitative information extracted from corresponding SHG images can be introduced in previously developed mechanical models to better characterize the behavior of arteries (Holzapfel et al. 2000;Morin et al. 2021). In the literature, researchers focused on three types of information that can be extracted from SHG images of collagen fibers: its geometry, its composition, and its morphology. However, they dealt with different types of input data (thus, the output data were different) at different scales of measure. In the present section, we will exhibit how that information was extracted through the literature. We summarize it in Table 4.

Geometric information (orientation, waviness)
A strong attention in the biomedical community is paid to geometric information of the collagen fibers. Changes in their geometric characteristics when they are under a mechanical load can actually be seen with a naked eye on SHG image; hence, the will to quantify it. Besides, their arrangement has a strong impact on the tissue's biomechanics.

Scale of measure
Orientation, waviness, and curvature are the important geometric information about collagen fibers. Orientation is usually calculated globally, but sometimes researchers focus on specific regions in an image and therefore on the local directions. On the other hand, waviness and curvature are determined locally.

Local characterization
The study of collagen fibers in biological tissues showed that those fibers are crimped and undulated. Thus, it is important to characterize their shapes. For this matter, several techniques have been proposed.
For the estimation of the fiber waviness, one needs to start by extracting the fibers. Sugita and Matsumoto (2017) determined the centers of the fibers as the pixels with a local maximum intensity. Then, they computed the length of the fiber as the distance between all the centers of a same fiber.
The CT-FIRE algorithm (Bredfeldt et al. 2014) is one of the techniques used to improve the images by extracting the fibers. The developers used also their algorithm to extract the collagen fibers and then estimated the waviness. Mathematical function -It fits well the orientation distribution profile of collagen fibers -It can be generalized to higher dimensions with few parameters -It assumes that the fibers follow one direction CT-FIRE was also used in Best et al. (2019) to extract the collagen fibers in renal cell carcinoma and by Zhou et al. (2017) in gastric cancer in order to characterize their organization and their straightness. It is also possible to segment the SHG images and extract the collagen fibers using other methods such as the skeletonization. Koch et al. (2014) proposed a new approach based on the application of several filters before segmenting the images. They used sequentially a CLAHE, a histogram adjustment, and a Frangi filter to reduce the noise and enhance the fibrous information. Then, a threshold was applied to recover a binary image where the fibers are well defined. Finally, they applied mathematical morphology operators to retrieve the fiber skeleton.
Techniques which were not initially developed for quantifying collagen in SHG images were also used. The most known one is the NeuronJ plugin of ImageJ software (Meijering et al. 2004). This plugin was designed to characterize neurons which have a linear shape. NeuronJ was used for tracing the fibers and analyzing their waviness (Zyablitskaya et al. 2017;Chow et al. 2014;Zeinali-Davarani et al. 2013). Besides, a 3D implementation of this technique was proposed and tested on SHG images. For example, to determine the fiber arc length, Hill et al. (2012) proceeded to a reconstruction of the SHG stack using a fast-marching algorithm to trace the fibers.
Once an accurate extraction of the collagen fibers is reached, it is possible to compute the waviness as a ratio between the Euclidean distance between the starting and ending points of a fiber and its actual length (Ayyalasomayajula et al. 2019;Hill et al. 2012;Koch et al. 2014). The estimation of those distances is done manually using ImageJ (National Institutes of Health, Bethesda, MD, USA) or Imaris (Bitplane, CT, USA).
The waviness in the 3D space was also investigated by Luo et al. (2017). They proceeded to a 3D skeletonization based on the fast-marching algorithm. The waviness computation is similar to what has been explained before, except that the considered points have 3D coordinates.
Regarding the local orientation, some interesting techniques were tested on collagen gels and showed their efficiency. One can cite the work of Bayan et al. (2009) where they used the Hough transform on different small partitions of the SHG image to determine the dominant local orientation of the considered fiber. The size of the partitions is chosen such as they are likely to contain a linear fiber. The SHG images were pre-processed to delete the noise through an adaptive thresholding and the application of an erosion and a dilation if needed.
It is also possible to evaluate orientations after fiber extraction. In Koch et al. (2014), the authors used the segmented skeleton to calculate the local orientation as the angle of the tangent line between the first and last points in a considered segment. Some other researchers used the FFT to evaluate the local orientation (Sivaguru et al. 2010;Rao et al. 2009;Ambekar et al. 2012a;Lau et al. 2012). For example, Rao et al. (2009) focused on the preferred orientation and the maximum spatial frequency of some regions in the SHG images. To determine those metrics, they computed the 2D FFT of the considered regions. The FFT gives the perpendicular angle to the preferred direction. To have a better quantitative approximation, one can fit the probability distributions of fiber orientations using one Gaussian function (Sugita and Matsumoto 2017). It is also possible to apply a 3D FFT on the entire stack to evaluate the fiber preferred direction in the space (Lau et al. 2012). However, the poor resolution of the SHG images in the third dimension may have a bad impact on the result of the 3D FFT to estimate the fiber directions in space.
Wavelet transforms were also used for direction estimation (Tilbury et al. 2014). The properties of the wavelet transform allow capturing small details and thus estimating correctly the orientation of the fibers. For this matter, the local coefficients of the wavelet transform were calculated and then clustered using K-nearest neighbors (K-NN) (Altman 1992) and principal component analysis (PCA) (Pearson 1901).
Image gradient is an efficient method to estimate orientations. This technique was initially developed by Chaudhuri et al. (1993). It consists in computing the gradient of the image to detect its edges and then to keep only the most relevant direction. The proposed method is similar to the Hough transform. It was later applied to biological tissues (Karlon et al. 1998) and to SHG images in particular (Hill et al. 2012;Phillippi et al. 2014;Kabir et al. 2013;Sun et al. 2015). In Kabir et al. (2013), the authors focused on a ROI from initial SHG image where the fibers have a pronounced dominant direction and calculated its 2D gradient to estimate the fiber orientation. A powerful ImageJ plugin that has proven its efficiency on biological images is Ori-entationJ. It is based on computing the image's gradient and its related weighted 2D structure tensors at each pixel. Cavinato et al. (2017) used this plugin to extract the orientation distribution histogram. Gaussian functions were then fitted to the histogram in order to quantify the dominant fiber directions. In Avila and Bueno (2015), the authors also used it on the image structure tensor.
Even though most of the proposed methods that have been used to quantify collagen fiber orientation were performed in 2D, some researcher such as Liu et al. (2018) took into consideration the collagen fiber distribution in the 3D space. They used the 3D directional variance algorithm to identify each pixel orientation and then estimate the entire fiber orientation.
More recently with the emergence of deep learning algorithms, some authors applied this technique to estimate local orientations of collagen fibers. For example, in Schmarje et al. (2019), a comparison of different 2D and 3D methods aiming at estimating local orientations was proposed. Besides, the authors introduced a new modality to transfer 2D weights to 3D weight in different-network architecture to perform a segmentation of some images with respect to local orientations.
Global characterization Most of the scientific contributions aiming at extracting quantitative information from SHG images of collagen fibers in biological tissues focused on the fiber orientation.
It is possible to determine fiber orientation using the FFT. It is the most used technique for this matter (Ayyalasomayajula et al. 2019;Bueno et al. 2013;Chiu 2010;Chow et al. 2014;Forouhesh Tehrani et al. 2021;Lau et al. 2012;Lee et al. 2019;Pijanka et al. 2019;Robinson et al. 2016;Sivaguru et al. 2010;Wu et al. 2011). This approach is also called FT-SHG imaging (Ambekar et al. 2012a). In Lee et al. (2019), the authors used the FFT on the entire stack and then performed a segmentation on the transformed images to only keep the dominant fibers' direction. Once the segmentation achieved, it is possible to recover the angles' distribution that corresponds to each image. It is then possible to evaluate the variation of the angles while going deeper in the stack. Germann et al. (2018) used the same methodology as Bueno et al. (2013), based on some pre-processing (noise reduction and edge sharpening) and a FFT to extract the orientation of collagen fibers in SHG corneal images.
Usually, the use of the FFT is sufficient to estimate the fiber directions, but sometimes it is useful to make the procedure more automated. For example, Ayyalasomayajula et al. (2019) extracted the distribution using a finite mixture of Von Mises distribution to fit the orientation distribution extracted from the FFT in order to determine the global mean orientation. Others, such as Schriefl et al. (2013) and Polzer et al. (2013), used a classical Von Mises distribution for the same purpose. It is also possible to use a Gaussian function for the fitting such as in Ambekar et al. (2012b). In some papers (Brisson et al. 2015;Kroger et al. 2021;Tang et al. 2014;Wu et al. 2011), the focus was oriented toward the result of the FFT where an ellipse was superposed. The major axis of this ellipse corresponds to the orthogonal of the dominant direction if the ratio between the major and the minor axes is high. Otherwise, there is no preferable orientation. Besides, it may be useful to use the radon transform on the 2D FFT of the SHG images (Mclean 2015;Mega et al. 2012) since, unlike the FFT, it provides quantitative information for each discrete angle. It is also common to use wedge filters after the FFT and then fit the orientation distribution with a Von Mises distribution to better estimate the orientation (Polzer et al. 2013;Schriefl et al. 2013;Niestrawska et al. 2016). In Zeitoune et al. (2017) the author applied an FFT on the images. Then, they improved the result of the transform by smoothing and enhancing it.
Directional filters were also used to determine the local orientation of the collagen fibers. Wen et al. (2014) proposed an approach based on those filters with different scales to determine the collagen fiber orientation in ovarian cancer. They extracted a histogram of the frequency of occurrence of individual patterns in an image. A nearest neighbor classification was then performed on the extracted histograms to distinguish between human normal and high-grade malignant ovarian tissues.
Some local techniques such as texture analysis have been used to quantify and describe the main fiber orientation. They showed their efficiency and they may be also more precise than the classic FFT. In fact, Hu et al. (2012) proposed a new approach for texture analysis based on orientation-dependent gray level co-occurrence matrix. They used their algorithm on ex vivo rat tendons to study the dominant collagen fiber direction. For this matter, they focused on the correlation feature of the GLCM.

Input data nature (2D/projected 3D/3D)
The determination of the orientation and the waviness of collagen fibers can be done using different types of input. Multi-photon microscopes allow going deeper in the tissue, and one can recover 3D stacks of images. However, most of the proposed techniques in the literature were limited to the image plane.
Generally, 2D images are used. For example, Zyablitskaya et al. (2017) used 2D SHG image of rabbit sclera to estimate the waviness of the collagen fibers. In addition, to assess the accuracy of their measurements, they calculated the average value on 10 SHG images. Ayyalasomayajula et al. (2019) used 2D images but limited their study to 10 slices of the stack. Then, the global orientation was set as the average of the computed orientations.
However, for the computation of the waviness of the collagen fibers, some papers processed 3D images such as in Hill et al. (2012), so they were able to characterize this metric in 3D in arterial tissues. In this paper, the waviness was computed from a 3D reconstruction of the SHG images by tracing the fibers using a 2D marching algorithm. This is possible because the waviness estimation is based on coordinates, which can be deduced from 3D SHG images. Meanwhile, an accurate 3D reconstruction of the SHG image may be hard to get because of the poor data resolution in the third dimension.
Regarding the orientation measurements, Hill et al. (2012) and Cavinato et al. (2020) used a 2D superimposed projection of the 3D stack of SHG images. Phillippi et al. (2014) succeeded in evaluating both collagen and elastin fibers in the aorta using superimposed 2D image stacks. The dominant orientation from a projection of all the SHG stack images can be extracted to recover a 2D image that contains all information from the entire stack (Hristu et al. 2018). However, it is more common to use 2D images to estimate the orientation (Bueno et al. 2013;Kabir et al. 2013) and look at its evolution with respect to the stack depth (Lee et al. 2019).
Some studies showed that the collagen fiber orientation in the axial-radial direction is negligible (Humphrey and Holzapfel 2012;Wagenseil and Mecham 2009). However, in Lau et al. (2012), the authors proposed a 3D FFT approach to evaluate the fibers preferred orientation in 3D stacks of SHG images. SHG stacks were also used to determine the waviness of the fibers such as in Luo et al. (2017). Bivariate Von Mises distribution was also used on 3D stacks of collagen in the aorta to fit the in-plane and out-of-plane collagen fiber orientations (Niestrawska et al. 2016).

Output data nature
The outputs of all the methods cited above can be divided into two types: a single value or an orientation distribution (i.e., a list). For example, the use of the FFT followed by ellipse fitting (Brisson et al. 2015;Tang et al. 2014;Wu et al. 2011) gives one value corresponding to the dominant orientation in the considered stack or ROI. Single-orientation values can also be extracted using texture analysis, which is applied locally (Hu et al. 2012). It is also possible to extract an orientation distribution histogram from the FFT by the application of a radon transform (Mclean 2015;Mega et al. 2012). It is also possible to use a Von Mises distribution and fit it to the orientation distribution obtained by the FFT (Ayyalasomayajula et al. 2019;Niestrawska et al. 2016;Polzer et al. 2013;Schriefl et al. 2013). Histogram of the frequency of occurrence is another representation of the orientation ).

Composition information (density)
Fiber density estimation is important for collagen characterization. In the literature, there are two ways to define the density: the volume occupied by the fibers in the stack (i.e., volume fraction) or the number of fibers in a considered region.

Scale of measure
The scale of measure depends on what has to be quantified. For volume fraction estimation, the procedure is global and applied to the entire stack. It may be interesting for some applications to focus on a ROI in the stack and calculate its volume fraction (for example, to characterize the evolution of tumors density). The same reasoning is applicable to calculate the density as the number of fibers in the entire stack or in a ROI. However, for an accurate estimation of the density, it is important to choose a ROI that covers up to 10 times the collagen fiber diameter.
In order to evaluate the fiber density, it is mandatory to enhance the SHG image by improving the signal-to-noise ratio. For this matter, it is important to filter the image and to recover an accurate representation of the fiber network through segmentation (Hompland et al. 2008) or fiber extraction (Wegner et al. 2017). Gade et al. (2019) performed a segmentation on the SHG stack using the Otsu thresholding. Then, the authors computed the area of segmented pixels in every slice and sum up the segmented area across the volume to calculate total areal density in the image stack. The same procedure was followed by Balu et al. (2014) and Tjin et al. (2014) where they performed a segmentation using ImageJ and then computed the collagen density as the sum of the pixels that have intensity values greater than a certain threshold.
It is sometimes interesting to proceed to a complete image enhancement step because of the diminution of the pixel intensity when we go deeper in the stack. For this matter, it is useful to apply a CLAHE on the SHG images. Cai et al. (2014) enhanced the dermal layer of human skin SHG images using the CLAHE algorithm. Then, they applied the Frangi filter and a segmentation using Otsu's tresholding in order to capture a representation of both the fibers and the holes in the images.
CT-FIRE is another algorithm used to extract the collagen fibers (Best et al. 2019;Wegner et al. 2017;Zhou et al. 2017). For example, in Best et al. (2019), the authors extracted the collagen in renal cell carcinoma in order to evaluate the density of the collagen in low-and high-grade tumors. The density can be calculated as the number of pixels corresponding to the fiber network with respect to the image or to the entire stack.
Second-order statistics, in general, and the grey level cooccurrence matrix, in particular, have been used to estimate the density of collagen fibers. In Kroger et al. (2021), the authors used the GLCM and especially the homogeneity parameter to determine the density of features in an image.
Some papers focused on the estimation of the ratio of both collagen and elastin fibers in SHG stacks (Abraham and Hogg 2010;Lin et al. 2005;Koehler et al. 2006). In Abraham and Hogg (2010), the authors started by filtering the images to reduce the noise. Then, they segmented the images and estimated the volume fraction of the fiber network as the sum of all the pixels belonging to the segmented region.

Input data nature (2D/projected 3D/3D)
The computation of the fiber density (also referred to as the volume fraction) requires the entire stack. The evaluation of the density can be done in 2D (i.e., slice per slice) or directly on the entire stack. It depends on how the segmentation is performed (in 2D or 3D). Cai et al. (2014) focused on 2D virtual biopsy images and not stacks. Therefore, they tested their approach only on single 2D images. In Zhou et al. (2017), the authors focused on single SHG images and, thus, calculated the collagen density in the 2D plane.
In general, SHG images are segmented separately. For example, in Gade et al. (2019) and Balu et al. (2014), the authors segmented the SHG images using a thresholding technique. Then, they calculated the number of white pixels in every image and summed them up across the volume to calculate total 3D density.
A global overview of the stack gives a more accurate estimation of the collagen density. That is why some researchers such as Abraham and Hogg (2010) implemented their method on the entire SHG stack.

Morphologic information (fiber's size)
In addition to geometric and composition information, it is necessary to know the fiber morphology in order to have a complete picture of the considered microstructure. For this matter, the intersections between the collagen fibers and the fibers size have been investigated in the literature.

Scale of measure
The study of the intersection between collagen fibers and even the estimation of their size is done locally because they are specific to a fiber (size) or a region (intersection). In order to be able to extract that information from the SHG images, it is important to enhance them.
For example, Koch et al. (2014) used segmented SHG images to estimate the fiber diameter. They performed some mathematical morphology operators (erosion and dilation) to obtain a 1-pixel-thick fiber skeleton. This skeleton was later used to calculate the fiber radii from the initial segmented image.
In some cases, depending on the application, only the characterization of the evolution of the morphology is needed. For this matter, texture analysis is used. In Wu et al. (2016), the authors used this technique to study the impact of aging on the skin microstructure. They computed the contrast, correlation, and entropy from the GLCM of the image and analyzed them to characterize the fiber structure and morphology. The contrast was computed to assess the presence of a fine structure of collagen fibrils. Wu et al. (2016) characterized how the collagen matrix is distinct from its surrounding and if there is loss in collagen through time using the computation of correlation. This can be generalized to distinguish between the collagen fibers and, thus, estimate their diameter (Cicchi et al. 2009). It is also possible to deduce if there are linear fibers and a fine structure through the computation of the entropy (Wu et al. 2016).
Some researchers showed interest in evaluating the fiber length. In Sugita and Matsumoto (2017), the authors extracted the centers of each fiber assuming that they correspond to a maximum intensity value and then estimated the fiber length as the sum of the distances between their centers. Fiber length can be evaluated manually from the SHG images of collagen gels after segmentation and using ImageJ drawing tool (Ajeti et al. 2011).
Moreover, collagen fibers that are extracted using the CT-FIRE algorithm can be used to extract manually the length and the fiber diameter (Drifka et al. 2016;Rosen et al. 2020;Wegner et al 2017;Zhou et al. 2017). In Rosen et al. (2020), the collagen fibers in every SHG image of feline mammary adenocarcinoma were identified by the mean of the CT-FIRE algorithm. Once the fiber extraction is achieved, each fiber was analyzed and its length and width were extracted in addition to the percentage of straight fibers.
Some out-of-the-box techniques were used. For instance, Robinson et al. (2016) estimated the collagen fiber thickness using the BoneJ plugin (Doube et al. 2010) for ImageJ which was initially developed to measure bone geometry. This algorithm gives the thickness of a considered fiber.

Input data nature (2D/projected 3D/3D)
For texture analysis, the application is usually performed on 2D images. Indeed, Wu et al. (2016) were only interested in investigating some layers of the dermis with the strongest collagen intensity. In Cicchi et al. (2009), the investigation was also limited to 2D SHG images of human dermis.
In the case of a skeletonization, the authors of Koch et al. (2014) used two 2D images (one of the fiber skeleton and one of the enhanced image) to determine the fiber radii. In addition, to estimate the fiber length, they used the skeleton of the fibers. Sugita and Matsumoto (2017) focused also on the fiber length and used 2D SHG images since fiber centers were determined in the 2D plane. It is also possible to extract the fiber network using the CT-FIRE algorithm and determine the fiber length and diameter (Drifka et al. 2016;Zhou et al. 2017). The papers that considered segmentation of the SHG images focused on each image individually and did not apply the segmentation to the entire stack (Ajeti et al. 2011).

Comparison of some methods
In the following, some of the collagen quantification methods that were previously described will be tested on a case study SHG image. The choice of those methods is based on their efficiency and how often they are used. Therefore, the tested methods are the FFT, the gradient through the Orien-tationJ plugin, and the CT-FIRE. The needed pre-processing is described. We performed all the cited methods on an SHG image of the adventitia layer of a human aorta , Fig 2.a. In the following, the estimated angles are expressed with respect to the horizontal direction.

FFT
The FFT was applied on a SHG image of collagen fibers of human aorta. We first used the FFT on the entire image to extract the dominant orientation in the considered image. Then, for comparison purposes, we focused on a ROI of the image where the fibers are aligned with respect to one direction.
The FFT of the SHG image was computed and a thresholding operation was applied on the FFT result. For this purpose, we used ImageJ. Then, an ellipse was fitted on the resulting image to recover the main direction of the collagen fibers. It corresponds to the direction perpendicular to the angle of the major axis of the ellipse. Figure 2.b exhibits the result of the FFT of the considered SHG image after thresholding. In our case study, the dominant orientation of the collagen fibers corresponds to 41.954°.
We applied the FFT on a ROI from the initial image. The ROI was chosen as an area containing only one fiber. The ROI was smoothed using a median filter, Fig. 3.a. Then, the FFT (Fig. 3.b) and the power spectrum (Fig. 3.c) of the ROI were calculated. As can be seen from both the FFT and the power spectrum representations, the result emphasizes one angle that corresponds to the longest portion of the fiber. Our calculations give an angle around 61° as the dominant orientation.
One can clearly see that in the presence of a distinct orientation, the result of the FFT gives a good estimation of that orientation.  (Hill et al. 2012;Phillippi et al. 2014;Kabir et al. 2013;Cavinato et al. 2017

Gradient
To calculate the gradient of the SHG image presented in Fig. 4.a, we used the OrientationJ plugin on ImageJ. First, we tried a global approach to detect the dominant direction of all the collagen fibers. Figure 4.b shows the result of orientation distribution. One can see that the preferred orientation of the fibers is around 42.5°.
To limit the non-relevant calculation due to noise, we applied a median filter on the initial SHG image to smooth it and to reduce the artifacts inside the fibers. Figure 5 shows the filtered image and the new orientation distribution. The orientation distribution presents one major peak located around 43°.
The gradient method was also applied on the ROI previously considered. The results can be seen in Fig. 6. The orientation distribution shows a dominant peak at the angle 66.5°.

CT-FIRE
The CT-FIRE algorithm was first applied to the initial SHG image without any pre-processing. The result is exhibited in Fig. 7. Because of the poor quality of the considered image, the CT-FIRE could not provide a good extraction of the collagen fibers.
As the fiber extraction was not good, the estimation of the global orientation (63.4°) is far from what has been computed using the FFT and the gradient (around 42°). The same steps were applied on the smoothed version of the SHG image and the results were similar. Indeed, the orientation was estimated to be equal to 61.3°.
We then applied the CT-FIRE algorithm on the same ROI. The algorithm was able to extract the skeleton of the fiber, Fig. 8. It estimated the fiber orientation to be around 70°.

Comparison
The experiments conducted previously showed that the FFT is an appropriate tool to estimate the main orientation of the collagen fibers if the fibers are well organized. Otherwise, no distinct information can be retrieved from the FFT. Besides, the result is highly sensitive to the thresholding algorithm used. On our chosen SHG image, the Otsu thresholding did not give an accurate orientation estimation. Moreover, for noisy images, the result of the FFT may be noisy too and, thus, non-exploitable.
The study of a ROI containing one fiber with two orientations and different intensities showed that the FFT focuses only on the most "visible" portion of the image. Our tests allowed us to extract only one angle from the FFT.
The global overview of the collagen fiber orientation using the gradient gives a result close to the one found using the FFT. This can be explained by the fact that the fibers in the considered image are not very crimped and undulated. Because the gradient is locally calculated, it is very sensitive to the shape of the collagen fibers. Therefore, it considers the fiber geometry. This method will always provide an estimation of a main orientation even if this one does not really exist. The result of this method is highly dependent on the quality of the image. Indeed, since it is a pixel-wise technique, it depends on the difference between the neighboringpixel intensities that can change while filtering the image.
When computed on a smoothed ROI with just one fiber, the gradient gives a good estimation of the orientation since it only involves information relevant to the considered fiber.
Regarding the CT-FIRE algorithm, its application on a bad-quality image without any pre-processing showed its limitation in extracting the collagen fibers. In addition, this technique looks like it computes the average of all the local orientations and not really the main global direction. The estimated orientation given by the CT-FIRE algorithm is close to the orientations given by the FFT and the gradient on a particular ROI of the image.
It is, however, interesting to use the CT-FIRE algorithm on a ROI. The experiments showed that, for small filtered ROI where there is only one fiber, the algorithm is able to extract correctly the position of the fiber's skeleton and thus estimate its orientation and its width and length.

Discussion
To our knowledge, reviews dealing with the quantitative analysis of collagen fibers from SHG images are not very common. The only one that we identified treated the topic from a method point of view. In fact, the main contribution of the present paper lies in the categorization of the image processing method with respect to the information that we want to extract (geometry, composition, or morphology). This structure makes it easier to biomedical researchers to find the most suitable method to the problem they are trying to solve. On another hand, the third part of the present review, which compared three of the most used techniques to estimate collagen fiber orientations, already gives the user an idea about how to use those methods and what to expect from them. Many image processing methods have been used in order to extract valuable information from collagen SHG images. However, the choice of the methods to choose depends deeply on the quality of the images and, thus, on the used microscope. In most of the cases, a pre-processing phase is necessary. In addition, the result of the pre-processing can affect, for example, the dimension estimation. In fact, if the image is not well filtered, blur can be mistaken to be part of a fiber.
It may be interesting to use machine learning algorithms to quantitatively analyze SHG images of collagen fiber. Some encouraging attempts can be found in the literature. For example, in Schmarje et al. (2019), the authors used convolutional neural networks (CNN) to segment SHG images of collagen fibers in order to quantify their local orientations.

Conclusion
In order to study the collagen fiber behavior in biological tissues, it is necessary to extract quantitative information to characterize them. To analyze those fibers, several techniques (including pre-processing and information selection) can be used, each one of them having advantages and drawbacks. The choice of the method to use is highly dependent on the information that need to be quantified. In this paper, we exhibited the most used techniques to quantify collagen fibers and we discussed the types of information that they allow to extract. As an illustration, we proposed a comparison of implementations of some of these methods to discuss their actual abilities to quantify collagen orientation. The choice of the method still depends on the images that need to be processed, their quality, and the error tolerance rate. A proper quantitative analysis of collagen fibers needs a combination of some of the techniques presented previously.
On the other hand, the quantitative analysis of collagen fibers in 3D is still not widely developed because of the limitations of the acquisition technique when going deep into the tissue and the poor imaging resolution in the third dimension. Further studies need to be oriented toward this issue especially because it is important to quantify the fiber network in the 3D space.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.