Hessian analysis for the delineation of amorphous anomalies in optical coherence tomography images of the aortic wall

: The aortic aneurysm is a disease originated mainly in the media layer of the aortic wall due to the occurrence of degraded areas of altered biological composition. These anomalous regions affect the structure and strength of the aorta artery, being their occurrence and extension proportional to the arterial vessel health. Optical Coherence Tomography (OCT) is applied to obtain cross-sectional images of the artery wall. The backscattering mechanisms in tissue make aorta images difficult to analyze due to noise and strong attenuation with penetration. The morphology of anomalies in pathological specimens is also diverse with amorphous shapes and varied dimensions, being these factors strongly related with tissue degradation and the aorta physiological condition. Hessian analysis of OCT images from aortic walls is used to assess the accurate delineation of these anomalous regions. A specific metric of the Hessian determinant is used to delineate degraded regions under blurry conditions and noise. A multiscale approach, based on an anisotropic Gaussian kernel filter, is applied to highlight and aggregate all the heterogeneity present in the


Introduction
Aortic aneurysm is a cardiovascular disease affecting the wall of the aorta artery. Degradation of the wall consists mainly of the loss of smooth muscle cells and elastin fibers, inducing weakening of the vessel wall [1]. The main manifestation of the pathology is a vessel diameter increase, continuously undermined as the periodic heart pumping softens its structure. When the artery is too weak to stand this heart pump rushing, it fissures producing a massive blood leakage. Before getting to this point, surgical intervention is needed to prevent patient's death. Surgical interventions consist on placement of an endovascular stent-graft or, when the pathology is severe, an external graft prosthesis by open surgery intervention. In the last case, the part of affected tissue is excised and replaced by a synthetic tube. The suture must be sewn to a healthy region of the remaining vessel or, otherwise, the aneurysm will reproduce or cause an aortic dissection.
The dimension of the areas of vessel wall affected by degradation is the real magnitude of degradation [1]. Hence, an accurate quantification of this magnitude is a hallmark of vessel health. To assess the actual condition of the artery wall, conventional staining techniques (such as Hematoxylin and Eosin or Verhoeff's Van Gieson) combined with common microscopy allow the visualization of the status at a cellular level. However, these techniques are time consuming and difficult to be integrated into the surgical workflow. Moreover, since the aortic tissue is thick and difficult to handle, the histological evaluation is limited to specific sampling points, once the tissue is excised, losing the evaluation of overall spatial status of the vessel wall.
Optical Coherence Tomography (OCT) is an imaging technology that allows the visualization of tissue samples in depth, imaging few millimeters within the sample. Resulting OCT images from tissue are equivalent to those obtained with common histology techniques, with the advantages of continuous and fast acquisition, easy incorporation into the surgical procedure but at fairly less resolution. Previous work [2] shows how the OCT approach correlates with the information provided by common histological techniques. Therefore, OCT could provide real-time guidance information about the vessel state under intraoperative conditions if an accurate and automated delineation of degradation is developed.
For delineation purposes, Hessian analysis has been previously proposed to detect vascular vessels, with tubular shape, in two [3,4] and three dimensional [5] images obtained by Magnetic Resonance Angiogram (MRA). Further modifications of this method were developed to detect nodules and vessels, based on the detection of predefined lines and bloblike [6,7] or curvilinear structures [8]. Multiscale analysis [3] has been also used to analyze OCT images for the study and quantification of vessel shape in the case of lymphatic vessels [9], angiography in the retina [10] and in extremities [11] and for segmentation purposes in the case of pig airway wall [12].
In contrast with previous experiences of Hessian analysis, in the present study of degradation of the aortic wall, the main issues are the lack of predefined shape and size of the pathological regions and the presence of high levels of noise due to reflections in small inhomogeneities along the wall thickness. In addition, attenuation of the aortic tissue strongly reduces sensitivity within depth, reducing the dynamic range of deeper areas (in contrast with the more superficial ones) and blurring anomalies boundaries, making anomaly segmentation a challenging process [3,13,14]. An alternative Hessian analysis is proposed to delineate degradation areas of amorphous shape, different sizes and intended to work under noise and variable dynamic range conditions, as it occurs in OCT images of the aorta wall. A new delineation accuracy estimator metric is also proposed to evaluate and optimize the anomaly detection performance avoiding a subjective evaluation. Finally, the detected anomalies have been used to quantify degradation in the aorta and this degradation has been compared against its histological quantification.

Aorta samples and protocol
In this viability study, human ex-vivo samples were obtained from the thoracic ascending aorta (TAA) region. Specimens come from two groups for comparison: healthy donors (11 patients), in compliance with the standard heart donor criteria, and diseased aneurysm samples (21 patients) from aneurysm surgical open repair interventions. Specimens were conserved in Phosphate Buffered Saline (PBS) and measured with the OCT system within 12 hours from excision. Afterwards, they were conserved in a diluted formaldehyde solution upon histological analysis.
OCT measurements were taken at several regions on each specimen and 5 OCT images, i.e. B-scans, were taken on each regions (356 B-scans for the whole specimen set). The number of regions depends on the size of the specimen, being a number between 1 and 4, with dimensions of 10x10mm. During histology, several cuts were applied to the regions, considering only the quantification of the region presenting the highest degradation.

Histological analysis
Five different stains and procedures were applied to each sample: Hematoxylin and Eosin, Verhoeff's Van Gieson, Alcian Blue at pH 2.5 and Alpha-smooth muscle actin [2]. The dimensions and severity of degradation in the aorta wall were measured by expert pathologists, providing the degradation state of the affected regions, quantified with a degradation histology score [1,15], grading from null degradation (score 0) to severe degradation (score 15).

Optical coherence tomography setup
All samples were analyzed with the Thorlabs OCS1300SS swept source system. The working wavelength is 1325nm ± 50nm bandwidth. The sensitivity of the system is up to 100dB. Cross-sectional B-scan images have dimensions of 10mm width by 3mm depth in air, sampled with 1024 pixels width by 512 pixels depth. Intensity level was digitalized with 8 bits.
The imaging probe was focused on the intima surface of samples to obtain information of the first layers seen from the lumen, comprehending the whole intima layer and a section of the media layer. According to OCT performance, light is backscattered when a refractive index change appears within the sample. The received signal is then proportional to the refractive index variation. The aortic wall tissue is turbid and thick, producing numerous small inhomogeneities and attenuation along penetration, making OCT images noisy and blurry. These effects reduce the original system dynamic range to approximately 40 to 60dB.

Amorphous spots identification
For the detection of specific anomalies, it is necessary to have previous knowledge about size, shape and intensity of the region surrounding each pixel and the global intensity of the image. The procedure to detect amorphous regions is divided into different sequential steps (Fig. 1). • Air-tissue interface segmentation: the air region between the optical probe and tissue must be removed prior to the analysis. In this work, this region is eliminated by applying a strong median filter (30 pixels vertical by 60 pixels horizontal) to the OCT image and posterior Otsu's thresholding algorithm [16] providing the air-tissue interface. The filtered image is discarded and only the air-tissue interface is preserved. The raw image below the air-tissue interface is not affected by this filtering process.
• Hessian matrix: the Hessian matrix of an image is defined as follows for every pixel f(x,y) in the image: As described in Eq. (1), for every pixel f(x,y) in the image, the second order partial derivatives must be obtained. The chosen differentiation method is the convolution with the second order partial derivatives of a selected kernel, k(x,y). As convolution is a linear method, the procedure of kernel differentiation and following convolution with the desired image, is equivalent to applying a convolution (filtering) followed by differentiation of the filtered image resulting from the convolution. This is depicted in Eq. (2).
( ) The relation described in Eq. (2) is only true when functions f(x,y) and/or k(x,y) are differentiable, and can be applied for a second order differentiation, third order, etc as long as this condition is satisfied. The main advantage is that only the Gaussian function has to be differentiated, instead of differentiating every image.
As mentioned previously, the anomalies observed in OCT images of the aortic wall present amorphous shapes of different sizes and axial directions, being frequently horizontally elongated [2]. As there is not a relevant geometric function to model these amorphous shapes, the two dimensional Gaussian function in Eq. (3) has been assumed as the kernel for the derivation, k(x,y), also with smoothing and border enhancement purposes.
The second order derivatives used for convolution are presented in Eq. (4) and its graphical representation, dependent on the value of the kernel parameters, is shown in Fig. 2.
being σ x and σ y the standard deviation in the x and y dimensions respectively and µ x and µ y the Gaussian mean values. The kernel size in pixels is defined as 6 times the biggest σ x or σ y , in order not to truncate the kernel response. The decomposition of the Hessian matrix can be expressed according to the following expression: where the solutions are the eigenvalues λ 1 and λ 2 and the eigenvectors v 1 and v 2 that satisfy the equality condition, being λ 1 <λ 2 . The signs of the eigenvalues show the evolution of the backscattering intensity signal of adjacent pixels. This procedure is applied to each pixel on the image to be analyzed obtaining a pair of eigenvalues (λ 1 , λ 2 ) and eigenvectors (v 1 , v 2 ) for each pixel.
• Hessian analysis: the relation between λ 1 and λ 2 provides information about the shape of the region in the image [3]. For the detection of amorphous shapes, only λ 1 , which represents the direction of less variation on the backscattering signal, is needed. The sign of the eigenvalues is due to the intensity of the region. By selecting only negative eigenvalues, λ 1 <0, dark spots can be detected, which is the case study of this work. For every possible i combination of (σ x , σ y ), the obtained λ 1 is normalized, obtaining the parameter A i (σ x , σ y ): This parameter is computed for every pixel, accounting for intensity variations in its neighborhood.
• Anomaly delineation: when the Hessian analysis is applied pixel by pixel, its response is limited to the size of the applied Gaussian kernel (defined as the 6 times the biggest σ x or σ y ). For this reason, after the application of the Hessian analysis, Otsu's method [16] is applied to evaluate the intensity relevance according to the global image intensity histogram. Otsu's thresholding procedure is used to automatically delimit the anomalies according to the A i values in Eq. (6). This procedure finds the optimum threshold to binarize the image according to the histogram distribution of the processed eigenvalue in the whole image. As a result, a mask is generated with the delimitation of intensity anomalies detected in the image.
• Accuracy estimator metric: edge detection accuracy becomes difficult to evaluate without further information about the real anomaly extension. A gold standard segmentation procedure, such as manually delineation, is not feasible in terms of time consumption, automation, and is also loose due to faint anomaly edges present in the sample set images.
A preliminary evaluation of the attained delineation accuracy is done by visual inspection, but this method presents subjectivity. To supplement this method, a new accuracy estimator metric, Δ, is proposed in Eq. (7) to evaluate the best delineation performance.
Anomalies Enclosing After the delimitation of anomalies (Fig. 3), the mean intensity of detected anomalies ( Anomalies I , Fig. 3(b)) is computed and compared with the mean intensity of the area that surrounds and encloses the anomalies ( Enclosing I , Fig. 3(c)). The dimension of the enclosing area is defined by a width σ x (Fig. 3(c)).The difference of mean intensities of those areas is used as the accuracy estimator. When the delineation is accurate, the anomaly will be ideally dark and the surrounding ideally bright, presenting big differences in the mean intensity. Delineation errors will reduce this mean intensity difference. As real anomalies are not perfectly bright or dark, this metric is only viable for comparison of performance in the same image.
• Multiscale aggregate: as stated, the response of the differentiation and filter kernel is a function of the standard deviation (σ y and σ x , manually selected by the user according to an approximate expected anomaly size) and the center point (µ x and µ y , equal to zero for filtering purposes). The kernel width is defined by the standard deviation, which is related to the dimension of the anomaly on x and y axis. In addition to the dimension of the standard deviation, the ratio among σ y and σ x has been used to evaluate different shapes.
In the case of ratio = 1, the spots to be detected tend to a circumferential shape. When ratio>1, σ y >σ x and according to Eq. (4), a horizontally elongated derivative kernel function is applied. When ratio<1, σ y <σ x , a vertically elongated kernel is employed.
The delineation obtained in previous steps is strongly dependent on the values selected for σ y and σ x and the different ratio between these two values. Following a multiscale procedure, the delineation process is repeated for different σ y and σ x values, obtaining different anomaly quantification, A i (σ x , σ y ), in each case. The final result is the one that maximizes the anomaly quantification metric, A o : In order to produce the best combination of different A i (σ x , σ y ) values, the accuracy metric Δ defined in Eq. (7) is used for automated optimization of the parameters in the multiscale aggregate. It works as follows: A i (σ x , σ y ) values are introduced in Eq. (9) one by one; if the addition of a new A i (σ x , σ y ) provides a bigger Δ accuracy metric for A o , this A i (σ x , σ y ) is added to the A o aggregate; if the Δ metric for A o decreases, this A i (σ x , σ y ) is discarded.
• Degradation quantification: The area of anomalies is strictly linked to degradation in the sample [1]. A quantification score has been developed to assess degradation quantification. For each B-scan, degradation is computed as follows: where a k is the area of every individual anomaly in pixels and Na represents the number of anomaly areas in each B-scan. The exponential provides an approximation to what is visually seen: few big anomalies are worse than many small areas in terms of degradation [1]. The number of total pixels (512*1024) in an image is used to normalize the exponent. When the degradation of the whole set of B-scans for a single specimen is aggregated, the OCT score for this specimen can be estimated: where Nb represents the number of B-scans for the sample. The term 1/Nb is included to normalize the OCT score, being independent of the number of B-scans of each specimen.

Results and discussion
The algorithm has been applied for the delineation on the whole set of aorta specimens described in section 2.1. For discussion purposes, only one of the worst cases of degraded aorta is displayed. Its corresponding OCT image is also affected by presence of high noise and tilt (Fig. 4) so, the validation of the method here assures its behavior under better conditions. The overall performance of the diagnostic procedure will be validated against the gold standard provided by histopathological analysis. The first steps of the analysis comprehend the air-tissue interface detection and the computation of the Hessian matrix of the image, resulting λ 1 and the consequent A i (σ x , σ y ) (Fig. 5). Different σ x and σ y will produce different kernel sizes, highlighting anomalies of different sizes and shapes. Fig. 5. Individual A i (σ x , σ y ) obtained for different σ x and σ y parameters. Pseudo-color represents normalized A i (σ x , σ y ) from its minimum value (blue) to the maximum (red) with 255 colors representation. Pink lines represent the air-tissue interface computed in the preprocessing stage. Application of Otsu's threshold method to these A i (σ x , σ y ) provides the anomaly delineation (black lines).
In a first analysis, the dependence on the kernel size is evaluated. Different σ x values have been applied to the selected image ( Fig. 4(b)), considering an isotropic kernel σ x = σ y , i.e. ratio = 1. When σ x is small, noisy areas are detected (Fig. 6(a)), whilst for large σ x the borders of the spots are lost and delineation becomes less accurate (Fig. 6(d)). In this case, intermediate σ x values provide a better delineation (Fig. 6(b), Fig. 6(c)). Both, visual inspection and accuracy metric Δ confirm this conclusion. The optimum result in this case, σ x = σ y = 7 pixels, is only valid for the OCT image under analysis. The optimum value can vary from one B-scan to another as a function of its intensity variations. When the ratio among σ x and σ y differs from 1 (Fig. 7), the response of the algorithm becomes biased towards the detection of vertical or horizontal regions. In this study, ratios below 1 remark vertical features in tissue depth ( Fig. 7(a)) that, in the case of the aortic wall result in wrong delimitation. On the other hand, ratios greater than 1 (Fig. 7(c)) can improve the identification of degraded areas as they tend to be horizontally elongated because the wall degradation related to the aneurysm disease is originated by the emergence of disorders in the horizontally oriented fiber structure of the media layer [2]. The tilt of the sample surface, and consequently the anomalies orientation, does not affect the accuracy performance for ratios around 1, but for large ratios (Fig. 7(d)) it results in a bad delineation. Fig. 7. Analysis of the delineation dependence on the kernel size orientation maintaining σ x = 7 always fixed and varying σ y as a function of the ratio value. Green lines represent the algorithm delineation output. Blue shaded regions represent manual reference delineation.
As seen in Fig. 6 and Fig. 7, the anomaly delineation performance of the Hessian analysis depends on its initial parameters. The determination of the best σ x and σ y is done by visual inspection of the delineation accuracy, numerically assisted by the accuracy metric in Δ (Table 1). In this case, intensity is computed from 0 to 255 levels, so the metric is expressed in this scale. According to the accuracy metric on Table 1, the best delineation, i.e. the one that achieves the biggest metric, is the case σ x = 7 and ratio = 1 (Δ = 32.4). The worst result is obtained for σ x = 12 and ratio = 3 (Δ = 16.9). The visual approach to the information exposed in Table 1 is shown in Fig. 8. The best simple implementation, i.e. without a multiscale approach, is shown in Fig. 8(a) and compared with the case where the accuracy metric Δ has been used for feeding the multiscale aggregate ( Fig. 8(d)). When Eq. (9) is applied to generate the group for the multiscale aggregate as described in Section 2.4, the best performance is achieved for a combination of σ x = 3, 5, 7 and ratio = 1 (Δ = 34.7), improving the case of a simple implementation σ x = 7 and ratio = 1. Fig. 8. Comparison of the best single application of the algorithm (top, σ x = 7, ratio = 1) and the best multiscale aggregate according to Δ (bottom, σ x = 3,5,7 and for each of them ratio = 1). On the first column (a, d), the original image with anomalies remarked in green and enclosing area in blue. The second row (b, e), present the segmentation of the anomalies with its mean intensity (I). The third column (c, f) represents the regions that surround the anomalies, required for the estimation of the accuracy metric, and its mean intensity (I).
Manual delineation is a laborious and not obvious task because it is affected by bias and unaccuracies, even in the case of trained users, that cannot neither perceive blurred nor small intensity variations. When automated delineation (Fig. 9(c)) is compared with the manual one ( Fig. 9(b)) extra anomalies appear. If both cases are compared with the original image, before delineation ( Fig. 9(a)), it can be seen how manual delineation misses small dark anomalies, mainly due to the fact that they present diffuse margins. The different A i (σ x , σ y ) included in the final A o delineation ( Fig. 9(d) to Fig. 9(f)), provide segmentation of different types of anomalies. The smallest kernel ( Fig. 9(d)) delineates small inhomogeneities, being imprecise in the case of big anomalies. On the other side, big kernels ( Fig. 9(f)) miss small anomalies. The aggregate combines all individual segmentations providing the best delineation according to the metric Δ. The relevance of small or big anomalies is studied later and validated against the degradation magnitude from the physiological condition of the samples.
The final delineation of anomalies in two different OCT images from two aorta specimens is shown in Fig. 10, applied in the case of severely degraded (upper row) and lowly degraded (lower row) specimens. OCT images ( Fig. 10(b), Fig. 10(d)) and conventional histology images ( Fig. 10(a), Fig. 10(c)) are processed for anomaly delineation. Histological images, the gold standard, evidence the real degradation of the inner layers of the aortic wall. In the case of severe degradation, anomalies in OCT images are correctly delimited ( Fig. 10(a), Fig.  10(b)), presenting a wide extension. In the case of low degradation (Fig. 10(c), Fig. 10(d)), the delimited areas are smaller, and related to small inhomogeneities, in the surface (in the case of OCT) or at a cellullar level (in the case of histology). After the anomaly delineation, the delimited areas are processed to quantify the degradation of the vessel wall. The calculation of degradation on every B-scan according to Eq. (10), is displayed in Fig. 11 in a boxplot representation. For every specimen, each B-scan produce a slightly different B-scan score, giving a degradation fluctuation.  Figure 11 shows how B-scan results differ within the same specimen providing a statistical variation and manifesting that vessel degradation is not a continuous spatial phenomena [1,2,[13][14][15]. In the case of boxplots presenting a small deviation, e.g. specimen 22, there is a strong statistic, giving close degradation results for all the B-scans in the same specimen. In the case of boxplots presenting higher deviations (specimens 1, 5, 11, 15 and 21), the quantification differs strongly among B-scans suggesting strong degradation in different areas, but not along the whole specimen. When these samples are closely analyzed, it is observed that they correspond to severely degraded samples, as stated by the pathologists. The case of sample 11 is an exception; there was damage only in a portion of one of the regions of the specimen. During the histological analysis, this small region was not seen by the pathologists, and the rest of the specimen presented good health conditions. According to (Eq. (11), the whole B-scan information is summarized into a global OCT score for each specimen. This quantification is shown in Fig. 12(b). The degradation score provided by the pathologist's procedure is shown in Fig. 12(a) for comparison. The gold standard histology score provided in Fig. 12(a), shows how the origin of samples (donor or aneurysm) does not imply null or severe degradation. The cases of sample 1, 5, 15 and 21 show the highest histology score, meaning they present the biggest degradation of all samples. Donor samples present scores from 0 to 5, as well as many of the aneurysm samples.When compared with the information obtained by the OCT score, results are comparable. Severe degradation samples present a greater OCT score than other samples, and the rest of samples provide lower and near to each other results.
In the case of histological analysis, only one histological section was processed for each specimen due to the difficult handling of excised aortic tissue. However, in the case of OCT analysis, it comprehends various B-scans (5 to 20 B-scans per specimen). As degradation is not homogenous along and aside the whole aortic vessel, the quantification of degradation will differ among B-scans [1,2] due to this uneven spatial co-registration. However, severe degradation specimens, where patient's risk is high, are highlighted both in the histological and OCT analysis. This makes OCT a promising technique for in-vivo evaluation of aortic vessel degradation. Tissue sampling could be even higher in OCT, reaching one B-scan each 10 µm at a rate of 1.6 B-scans per second. Therefore, the non-invasive and fast acquisition of OCT images, jointly with its high spatial resolution over the tissue surface, makes this imaging technology a promising alternative to generate degradation maps of the aortic vessel able to guide surgeons during their interventions.

Conclusion
The delineation of anomalous regions in OCT images of turbid media such as the aorta artery is a demanding task. In case of aneurysm disease, and to be able to guide surgeons during real-time interventions, this delineation is required to assess about the degradation status of the aorta. Diseased vessels present amorphous shapes along the wall thickness whose diffused borders pose a limitation for the delineation purposes. A solution, based on the Hessian matrix analysis, has been presented providing accurate delineation results. The proposed procedure depends on specific parameters of the Hessian analysis that have been automatically optimized using an accuracy estimator metric that avoids any heuristic approach. After the delineation, the detected degraded regions are classified according to their degradation dimensions. This method is applicable to other types of images with small modifications.
The method has been compared with histological grading techniques, providing complementary results but with higher spatial scanning resolution in the case of OCT. The quantification of degradation provided by OCT exhibits a statistical variation among regions of the same specimen that evidences the disperse spatial evolution of the physiological phenomena of vessel degradation. OCT measurements provide a different approach to the interpretation of results. Firstly, OCT can be used as a complementary method during histological analysis. Statistical variations can indicate whether there is a strong fluctuation in the sample degradation per regions, hence, a more extended histological analysis (more regions under analysis) should be performed. Secondly, OCT shows capability for detecting strong degradation produced in the vessel wall. This feature, and the fast acquisition rate, make OCT systems a feasible substitute/supplement of conventional histology techniques during surgical intervention, being transferable to multiple areas of interest.