Enhancement of morphological and vascular features in OCT images using a modified Bayesian residual transform

: A novel image processing algorithm based on a modified Bayesian residual transform (MBRT) was developed for the enhancement of morphological and vascular features in optical coherence tomography (OCT) and OCT angiography (OCTA) images. The MBRT algorithm decomposes the original OCT image into multiple residual images, where each image presents information at a unique scale. Scale selective residual adaptation is used subsequently to enhance morphological features of interest, such as blood vessels and tissue layers, and to suppress irrelevant image features such as noise and motion artefacts. The performance of the proposed MBRT algorithm was tested on a series of cross-sectional and enface OCT and OCTA images of retina and brain tissue that were acquired in-vivo . Results show that the MBRT reduces speckle noise and motion-related imaging artefacts locally, thus improving significantly the contrast and visibility of morphological features in the OCT and OCTA images. and OCTA images. Comparison of the MBRT method with Frangi filter and Gabor filter showed superior performance of the MBRT in terms of preserving complex vascular and morphological features in the OCT images.


Introduction
Optical Coherence Tomography (OCT) is a non-invasive, high-resolution medical imaging modality that can resolve morphological features in biological tissue as small as individual cells, at imaging depths in the order of 1mm below the tissue surface , estimate and compensate the image aberrations, therefore can enhance the visualization of features of interest. These approaches require good phase stability of the system and dense scanning pattern, so their transition to the clinic use systems is difficult and filed of view is limited in in-vivo imaging.
The Frangi filter [34] approach is widely applied to enhance vessel-like shapes in 2D and 3D OCT images by using second derivative information of the images. However, its performance is strongly dependent on the presence of morphological features of different sizes (scales) in the same image. For example, the effectiveness of a Frangi filter designed for one scale can be affected by presence of morphological features of other scales. Moreover, the Frangi filter is very susceptible to noise contamination, especially at smaller scales (for example, capillaries in the human retina that have cross-sections of ~10 µm to 20 µm). Furthermore, the Frangi filter typically generates images artefacts that resemble small vessels in the OCT and OCTA images, which can compromise the quantitative analysis of blood vessel networks in OCT and OCTA images. Gabor filter is also used to enhance OCTA images; however, this method may artificially increase the size of the vessel [35], and similar to Frangi filter, different kernel sizes can interfere with each other as lacking of prior knowledge.
Multi-scale decomposition has gained high popularity in the recent years in physiological signal processing. Wavelet decomposition [36] uses wavelet transform to decompose the physiological signal into different scales, where the noise component is usually in the finest scale and can be isolated during signal reconstruction. One limitation of this approach is that the decomposition procedure is based on deterministic functions, therefore, basis-related artefacts could be associated with different wavelet function selections. Alexander et al. [37] proposed a new approach from Bayesian perspective to decompose the physiological signal into different unique scales, and its implementation on noisy electrocardiographical (ECG) signal showed significant noise reduction while preserving the useful signal features. This method, called Bayesian residual transform (BRT), has the advantage of isolating residual information at each scale from other scales, and each scale is calculated based on its own statistical characteristics.
Here we present, for the first time, enhancement of morphological and vascular features in OCT and OCTA images based on a modified version of the BRT approach (MBRT). Similar to ECG signal in 1D, in the MBRT approach, the 2D OCT original images are decomposed into several residual images at different scales. Afterwards, different weights are assigned to the residual images by a scale-selective residual adaption. Morphological and vascular features of interest are emphasized by empirically assigning larger weights to them, and the final enhanced image is assembled by summing up all the weighted residual images. In this manuscript, we describe the MBRT approach in detail and present enhanced images that were acquired in-vivo from human and rats with different customized OCT systems and a commercial system, and therefore with different spatial resolutions, image acquisition rates and noise levels. The vascular images from customized systems were also processed with a Frangi filter and Gabor filter, and compared with the MBRT enhanced versions.

Bayesian residual transform
The full derivation of one dimensional Bayesian residual transform has been presented previously [37], and a modified version (MBRT), based on two dimensional images, is presented graphically in Fig. 1. An OCT or OCTA image contains structures with different scales, from finest scale such as speckle noise, to coarsest scale, such as major blood vessels and retinal layers with relative uniform scattering potentials. Therefore, let A denote coordinate (x,y) in a image, and ( ) I A can be decomposed to a series of residual images ( ( ) ( ) ), with incremental scale, each of which represents information at a specific scale: denotes the summation of residual images with scales ranged from k to N (k≤N). Therefore, the equation can be also modified to: Computing µ 1 can be considered as solving an inversing problem: In general, residual image with k th scale information can be expressed as: which can be solved by a non-parametric Nadaraya-Watson kernel regression approach [38,39], and by using a Gaussian kernel with a scale k α , the inversing problem can be computed as following: Now an image is decomposed into several residual images at different scales, and since each residual image is independent to other residual images, and the original image can be reconstructed by simply summing up all residual images, i.e. inverse BRT. Structural enhancement was achieved by scale-selective residual adaptation: scale coefficient i η was applied to the corresponding residual image ( ) µ i A in order to enhance/suppress the structural information at the unique scale of the image: and the structure enhanced image is computed as the sum of scaled residual images: This non-iterative image decomposition using BRT method was using MALTAB with digital imaging processing package, with Nadaraya-Watson kernel regression (Eq. (6) written in C + + and called by MATLAB Executable, working on an Intel Core i7-4930K Processor PC with 16 GB RAM. The proposed structure enhancement was compared with the most used structure enhancement algorithms (Frangi filter and Gabor filter) on the OCTA images acquired from customized OCT systems. The algorithm was also applied on MATLAB, and the Frangi filter kernel sizes were set to match the sizes of the vascular in the retina (2 × 2 to 6 × 6 pixels) and brain (3 × 3 to 15 × 15 pixels). Gabor filter used Gaussian kernels with 9 equally spaced orientations [0-2π], 0.1 spatial frequency, and kernel sizes 6 × 6 to 60 × 60 pixels for rat retinal vascular images and 6 × 6 to 100 × 100 pixels for rat brain images. The Gabor filter local optimization scheme followed the method proposed by Hendargo et al. [35].
In order to provide quantitative comparison of the proposed MBRT method with Frangi and Gabor filters, two indices were introduced: dice similarity coefficient (DSC) [40] and vascular connectivity. The ground-truth vascular plexus was manually segmented, and then compared with binarized images processed by all three methods. The binarization thresholds were then automatically decided using Otsu's method [41]. The DSC was defined as: where TP, FP, and FN are the true positive, the false positive, and the false negative pixels. The vascular connectivity was calculated similarly to the approach proposed by Jia et al. [42]: the binarized images were skeletonized, and the vascular connectivity was defined as the ratio of the connected pixel numbers to the total pixel numbers in the skeletonized images.

OCT and OCTA images
The performance of the proposed MBRT method was first applied to a test image composed of Newton's rings of different sizes. Then the method was applied to a series of OCT and OCTA images acquired in-vivo from the rat retina and brain, as well as, from the human cornea and retina. The rat retina and brain imaging studies were approved by the University of Waterloo Animal Research Ethics Committee and adhered to the ARVO statement for use of animals in ophthalmic and vision research. The rat retina and brain images were acquired with a SD-OCT system that operates in the 1060 nm spectral region and provides ~3.5 µm axial and ~5 µm lateral resolution in retinal tissue with ~100 dB SNR for 1.7 mW optical power of the imaging beam [43,44]. Cross-sectional retinal OCT images with distinctive retinal layers were acquired to examine the algorithm's behavior on images with layered structure. OCT images were acquired from the rat retina at the ONH and its vicinity (3.4 mm × 3.4 mm; 512 A-scans × 512 positions × 4 scans/position) OCTA images were generated from this data by a maximum projection in depth. In a similar way, OCTA images were acquired from the tumorous rat brain (5 mm x 5 mm; 1000 A-scans × 1000 positions × 5 scans/position). The brain tumor was induced by injection of human brain tumor cells through a small craniotomy 18 days before imaging. The human cornea imaging study was approved by the Research Ethics Committee at the University of Waterloo and adhered to the tents of the declaration of Helsinki. The human corneal OCT images were acquired with a fiberoptic, SD-OCT operating in the 800 nm spectral range, that provided ~0.95 µm axial and ~4 µm lateral resolution in corneal tissue at 34,000 A-scan/s image acquisition rate [45]. Cross-sectional images (1000 Ascan × 4096 pixel/Ascan) were acquired in-vivo from the corneas of healthy human subjects. The human retina vascular images were acquired with a commercial swept-source OCT (SS-OCT) system (PLEX Elite 9000, Zeiss, Germany). The scanning area on the macula is 3 mm × 3 mm and 6 mm × 6 mm. The human retinal images processed by MBRT were not compared with other filters as the extracted images from the commercial system were already processed by the built-in, proprietary algorithm(s).

Results and discussions
To evaluate the proposed algorithm on a 2D image with multiple scales, a Newtown's rings phantom image containing multiple concentric rings of graded scales was created (Fig. 2(A)). A zero-mean Gaussian noise (variance: 0.7) was added on the original image ( Fig. 2(B)). The noisy image was decomposed to 5 different scales from finest to coarsest (Fig. 2(C)-(G)) and then reconstructed using BRT algorithm. It is easy to observe that the concentric ring structures were well described from scale 2 to 5, where the scale contains mostly the background noise. The final structure enhanced image ( Fig. 1(H)) emphasized the scales from scale 2 to 5, and the visibility of the concentric rings were largely improved from the noisy image. Next, an example of the proposed method based on modified BRT were implemented on the OCTA images of rat's retina (Fig. 3). Max projection over certain depth in retinal tissue (0-35 µm) was utilized to generate en-face OCTA images on the rat's retinal GCL + NFL layer. This layer contains first-order retinal arterioles and venues (size ~30-50 µm) and their branching capillaries with variable sizes (5-30 µm). The original image (Fig. 3(A)) was decomposed into 5 different residual images with even scale interval (Fig. 3(B)-(G)), where structures of large arterioles and venues were extracted from scale 4-5, and structure of small capillaries were extracted from scale 2-4. Scale one contains mostly the statistical noise information, mixed with only a little capillary structure information. Therefore, when the final image is computed, the weights are set to [0,1,1,3,1] from scale 1-5, respectively. The contrasts of the original and reconstructed images were adjusted so that the total image values are equal. Processing time on this OCTA image with size 1000 × 1000 pixels is plotted in Fig.  2(H) as a function of decomposed residual image numbers. Linear increase from 1.2 seconds to 2.6 seconds is expected with decomposed image numbers from 4 to 11.  Figure 4 shows qualitative comparison of the original and processed maximum projection OCTA images from the OPL, using both Frangi and Gabor filters as well as the MBRT to enhance the retinal vasculature. This OPL has the highest capillary density among three inner retinal vascular layers (GCL + NFL, IPL and OPL) and the diameter of the capillaries is relatively uniform, therefore, it is a good sample to evaluate the performance of the MBRT algorithm on enhancing image features of approximately the same size. The original image and the images processed with the Frangi and Gabor filters and the MBRT are shown in Fig.  4(A)-(D) respectively. Regions of interest in those images marked with the yellow rectangles are magnified for quantitative comparison in Fig. 4(E)-(H). In the original image, although the general vascular distribution can still be visualized, the useful vascular structure has low CNR, contaminated by speckle noise and the shadow artefacts. Frangi filter can pick-up the vesselness of the vascular image and enhance its contrast, however, from magnified view (Fig. 4(F)), it is easy to observe that the algorithm is underestimating the boundary of the vessel, and continuity of the vessels are not preserved. Both vessel diameter bias and discontinuity associated with the Frangi filter will jeopardize further vessel quantification, e.g. diameter measurement and perfusion density calculation. The Gabor filter preserves and enhances the vascular features by maximizing the vascular orientation. However, as seen in the magnified view (Fig. 4(G)), its performance on the looped and bifurcated structure is relatively poor, as the local optimization is applied to a single direction. In Fig. 4(H), it is easy to observe that the background noise is well suppressed, while the sharpness of the vascular boundary is preserved. Again, in the processed image, the non-blood perfusion area is much clear and the capillaries were well preserved and even enhanced.  Figure 5 shows the performance of a Frangi filter, Gabor filter and the MBRT for enhancing the visibility of the tumorous vascular structure in the rat brain. Figure 5 . The avascular area at the center of the images (white arrow) is the region of a cyst, where the scattering potential is different from the surrounding tissue. In the originate image ( Fig. 5(A)), the vascular structure is not pronounced and contaminated by both speckle noise and motion artefact, displayed as vertical stripes in the max projected image. Frangi filtered image overestimates and underestimates the vascular structure. Specifically, with complex perfusion distribution (Fig. 5(E)), the Frangi filter artificially generated some vascular structure that are not present the original image. The Gabor filter generally enhanced the visualization and continuity of the vasculature. However, the noise level is more pronounced and therefore obscures the small blood vessels. Moreover, as both filter iterations do not take into consideration the a priori knowledge, the small scale structures are interfering with the large scale structures. MBRT enhanced morphological images (Fig. 5(D), (H), (L)) shown improved visibility and the continuity of the vasculature. The vessels' boundaries are clear, which enables more precise measurement of the vessels size and density. The boundary of the cyst is sharper, allowing unambiguous segmentation of the cyst area (blue curved line). Overall, the proposed MBRT algorithm could be helpful in studies investigating brain tumor growth by offering a more precise evaluation of the morphological and vascular changes in brain tissue associated with the tumor.  Table 1 shows the DSC and vascular connectivity of images processed by all three algorithms (N = 9). MBRT method has on average 28% and 4% improvement of DSC compared to the Frangi and Gabor filters. Manal segmented images have almost perfect vascular connectivity (0.99 ± 0.01). MBRT preserves best the vascular connectivity among all three methods, and is on average 21% and 2% more efficient in vascular connectivity compared to the by Frangi and Gabor filters, respectively. Moreover, the computational times of Frangi and Gabor filters are on average ~1.1x and ~3.5x longer than of the MBRT. Next, the MBRT algorithm was applied to retinal images acquired from a healthy subject with a commercial SS-OCT system. Figure 6(A), (C) are the extracted images of human macula area with sizes 3mm × 3mm and 6mm × 6mm, and structure enhanced images are shown in Fig. 6(B), (D). It is easy to observe that the central avascular area boundary gets sharper while the structure of the vascular is well preserved. Enhanced contrast of major blood vessels will make further analysis, e.g. vascular segmentation and skeletonization, easier and unambiguous. Fig. 6. Enhancement of the vasculature in OCT images of the human retina using MBRT, acquired with a commercial SS-OCT system. (A,C) Extracted original images of healthy human macula with sizes 3mm × 3mm and 6mm × 6mm. (B,D) MBRT enhanced images. Note that the original OCT images were already pre-processed with built-in filters into the commercial SS-OCT system.
Finally, the feasibility of the proposed MBRT algorithm was tested on a series of crosssectional morphological OCT images. Figure 7(A) shows a large view of the rat's retina acquired near the ONH. Distinctive small features, including capillaries and external limiting membrane (ELM), and large features, like retinal blood vessels and choroidal vessels can be characterized in the original image. However, the layer segmentation algorithm would fail because of speckle noise and discontinuity of some thin retinal layers. Figure 7(B) shows the reconstructed image with much cleaner retinal layer and choroidal vessel boundary, and the magnified view of the retinal blood vessel area enables the unambiguous segmentation of a small blood vessel (red arrows in Fig. C and D). Original image and the reconstructed image using MBRT approach that represent the layered retinal structure near the ONH, shown as the yellow line in (C). The magnified view of the blood vessel region shows the unambiguous recognition of a blood vessel after image reconstruction using MBRT approach. (D-E) original and structure enhanced image using modified MBRT. Layered structure is improved so that the layer segmentation can be easily implemented afterward in (F). It's note that the basal cell layer is much more prominent after the structure enhancement and can be confidently isolated from the rest of the epithelium layer.
Note that the image is processed in such way that the features with large scales are enhanced, therefore, the small features, such as the ELM are suppressed. The in-vivo human corneal image (Fig. 7(E)) was acquired with an OCT system that utilizes a supercontinuum laser, with a noise level significantly higher compared to the noise of the superluminence diode (SLD) light source used to acquire the images from the rat retina and brain ( Fig. 7(A)-(D) and Fig. 5 respectively). Furthermore, the scales-associated noise varies along the depth, and as such, a careful selection of the MBRT coefficients is required to isolate the relevant morphological information from the variable unwanted scales. Figures 7(E)-(G) show the original and reconstructed image of the human cornea, respectively, while Fig. 7(G) shows the segmented cornea layers based on reconstructed image. It is clear that the proposed MBRT enhanced image improves the visibility of the different layers and thus allowing for easier and more precise manual segmentation, especially in the case of the basal cell layer where distinctive scattering pattern can be differentiated from the rest of the epithelium.

Conclusion
We have developed a modified residual Bayesian transform for image decomposition to multiple scales and enhancement of morphological and vascular features in the images, achieved by using a scale-selective residual adaption. Images acquired in-vivo from human subjects and test animals with a number of OCT systems with different spatial resolution and image acquisition rates, were used to test the general feasibility of the proposed algorithm on OCT and OCTA images. Comparison of the MBRT method with Frangi filter and Gabor filter showed superior performance of the MBRT in terms of preserving complex vascular and morphological features in the OCT images.

Funding
Natural Sciences and Engineering Research Council of Canada (NSERC, grant # 312037 and # 03966), Canada Research Chair (grant # 229272) and University of Waterloo Research Incentive Fund.