Quality improvement of OCT angiograms with elliptical directional filtering

: We present a method of OCT angiography (OCTA) data filtering for noise suppression and improved visualization of the retinal vascular networks in en face projection images. In our approach, we use a set of filters applied in three orthogonal axes in the three-dimensional (3-D) data sets. Minimization of artifacts generated in B-scan–wise data processing is accomplished by filtering the cross-sections along the slow scanning axis. A-scans are de-noised by axial filtering. The core of the method is the application of directional filtering to the C-scans, i.e. one-pixel thick sections of the 3-D data set, perpendicular to the direction of the scanning OCT beam. The method uses a concept of structuring, directional kernels of shapes matching the geometry of the image features. We use rotating ellipses to find the most likely local orientation of the vessels and use the best matching ellipses for median filtering of the C-scans. We demonstrate our approach in the imaging of a normal human eye with laboratory-grade spectral-domain OCT setup. The “field performance” is demonstrated in imaging of diabetic retinopathy cases with a commercial OCT device. The absolute complex differences method is used for the generation of OCTA images from the data collected in the most noise-wise unfavorable OCTA scanning regime–two frame scanning. demonstrates that the elliptical directional filtering correctly differentiates between noisy vascular networks and noisy avascular areas in the OCTA data obtained in the imaging of normal eye.


Introduction
Optical coherence tomography angiography (OCTA) is a set of methods derived from the optical coherence tomography technique (OCT) to enable imaging of the vascular systems embedded in the living tissues. Since the first report of the OCTA imaging in the human eye [1], an abundance of OCTA methods was produced by research groups around the world [2][3][4][5]. They all relay on a simple observation that the flow of the blood can be used to identify the location of vessels within the static tissue. This concept is a special case of using signal variations as a contrast mechanism for imaging of features exhibiting specific dynamic processes, within stationary objects. However, in living tissues, motion occurs within the entire imaged volume at different temporal and spatial scales. The goal of motion contrast methods is therefore to identify the temporal and spatial characteristics of the motion of interest and adjust the experimental procedures and data processing methods to highlight the areas where the dynamic processes of interest occur. The OCT technique is well suited to detect the motion of blood in the circulatory systems of the eye due to appropriate imaging speeds and resolutions. However, it requires adjustment of the imaging protocols to ensure the detection of fast (tens of mm/s) and slow (fraction of mm/s) flows in the arteries, veins and capillaries of the eye fundus. It has also limitations connected to noise generated by the optical, mechanical and electronic components. In addition, various unwanted eye motions are easily detected with OCT, too. Thus, the simple concept of blood flow as a contrast mechanism in OCTA visualization of the circulatory systems of the eye requires solving several practical challenges.
Mathematical methods used for identification of the OCT signal variations caused by the blood flow typically relay on statistical analysis of the series of data sets acquired over time from the same location in the object. The most often used motion detection metrics are variance or correlation of the amplitude and/or phase of the complex OCT signal (methods calculating differences are special cases of the above mentioned approaches) . To improve the outcome of those methods, i.e. to reduce noise in the OCTA images and fragmentation of the images of vessels, many data sets are often acquired over time from the same location in the object enabling application of various averaging techniques [26,[29][30][31][32][33][34][35]. As a result, noise is reduced and vessels appear less fragmented in the OCTA images. However, the consequence of the acquisition of larger data sets is longer imaging time and increased number of image artifacts caused by the involuntary eye motions. Therefore, smaller data sets are typically used in connection with numerical filtration methods of OCTA images to circumvent motion artifacts or to enlarge the imaged areas of the object.
Numerical filtration is always applied to OCTA images. However, the details and effects of these procedures are rarely described in literature. Generation of en face OCTA projections from stacks of C-scans, by maximum intensity projection or calculation of the axial mean value, is by itself a filtration method of the 3-D OCTA data sets. Often, median filtering, Gaussian image smoothing or some other denoising method is also applied to 3-D OCTA data or angiograms (OCTA en face projections) with empirically adjusted filter window sizes. Another typically used procedure is scaling (in practice always enlarging) of the images. This involves interpolation algorithms which act as image filters. Aside from using these customary methods, publicly available in image processing software packages, researchers have been also developing methods accounting for features characteristic of the OCTA angiograms. Two approaches frequently appear in the literature: Gabor filtering [29] and Hessian filtering [36,37]. Recently, a method of shape-based classification of "vesselness" with anisotropic Gaussian filter was reported for contrast and connectivity improvement of vascular images in OCTA [38]. Many of those methods were adapted for OCT angiography from other medical imaging techniques [39][40][41][42].
In this paper we introduce a method of elliptical directional median filtering of the 3-D OCTA data sets. The goal is to use a set of image filters to reduce noise in the OCTA en face projections and improve the continuity of imaged vascular networks. The core filtering procedure uses a concept of structuring, directional kernels of shapes matching the geometry of the image features. Similar idea was presented in [43], where the authors used square and triangle kernels for digital image filtering. Application of median filter along the vessels was presented in magnetic resonance imaging in [40]. Here, we use rotating elliptical kernels to match the orientation of vessels in OCT angiograms and apply median filtering within the best locally matching ellipses. The effect is generation of OCT angiograms with improved vessels continuity and reduced noise.

OCT angiography setups and methods
To develop the elliptical directional filtering methods presented in this article we first used OCT data acquired with the Revo NX device (Optopol Technology Ltd., Zawiercie Poland [44]). This spectral-domain OCT setup operates at 830 nm center wavelength, with 5 µm axial and 18 µm transverse imaging resolutions. The imaging speed, measured as line rate of the camera used for the spectral interference signal detection, is 110 kHz. The OCT angiography method incorporated in the device: SOA (spectral domain OCT angiography), uses the absolute complex difference OCTA technique [26] to visualize the retinal vasculature. Details of the scan protocols used for imaging of a healthy subject and three diabetic retinopathy cases, performed in an ophthalmology clinic, are given in Table 1. Angiograms were automatically generated as maximum axial intensity projections of the 3-D OCTA data using segmentation of the anatomical retinal layers. Images of the "superficial vascular plexus" (en face projection of radial peripapillary capillary plexus, superficial vascular plexus, and intermediate capillary plexus) and the deep capillary plexus are presented.
To prove that the approach is universal and can be potentially used with any OCT device equipped with OCT angiography data acquisition and processing methods, we applied it to images obtained with a research-grade spectral-domain OCT setup described in details in our article by Rumiński et al. [26]. The imaging was performed in the macular area of a healthy eye. Here, pairs of B-scans were acquired from the same location in the eye fundus to enable direct implementation of the absolute complex difference OCTA method [26] according to the formula: where I OCTA is intensity of the OCTA image points calculated from corresponding A-scans in the B-scan pairs. A are complex amplitudes obtained via Fourier transformation of the wavenumber-linearized and dispersion-compensated interference spectra. t stands for time and z, x are spatial coordinates in depth and along fast scanning axis. M is the number of Bscan repeats at subsequent positions in the object. Prior to calculations of the OCTA signal intensity, shifting of the A-scans and bulk motion correction was performed [26]. Details of the scan protocol are given in Table 2. In parallel with the OCTA images, structural cross-sections were generated as mean images of the B-scan pairs. They were used for the correction of the axial B-scans displacements by correlation. The corrections were applied to the OCTA cross-sections. Finally, the position of the outer plexiform layer (OPL) was identified in the structural Bscans by peak detection and a surface approximating the shape of OPL was fitted using Zernike polynomials, as described elsewhere [31], see Fig. 1(a). This surface was used to empirically select the depth boundaries for generation of the OCTA en face projections, i.e. OCTA angiograms. Fig. 1. Three-dimensional data sets. a. Structural data set used for segmentation of the outer plexiform layer. The fitted surface was used for selection of depth positions of en face projections as indicated by the yellow and blue lines. b. OCTA data set indicating the orientation of three planes of interest used in our data processing. Green outlines -B-scans, XZ planes parallel to the fast scanning axis; blue -YZ planes, parallel to the slow scan axis; pink -C-scans, XY planes, perpendicular to the scanning beam. c. OCTA en face projections of the retinal vessels and the inner capillary plexus were generated as maximum intensity projection between the yellow lines. An example angiogram is schematically shown with yellow outline. The angiogram of the outer capillary plexus was generated between the blue lines.

Elliptical directional median filtering
In the design of our image filtering method we used the basic knowledge of OCT imaging specifics and the architecture (anatomy) of the retinal vascular system. Our assumptions were as follows: -retinal vasculature is a web of meandering, locally highly directional tubes, oriented nearly perpendicularly to the imaging beam, thus vessels should appear as meandering, locally highly directional lines in OCTA en face projections, -noise and artifacts generated in B-scan-wise data processing are not correlated among OCTA cross-sections, while images of vessels are correlated among neighboring Bscans, -depth resolution can be sacrificed within reasonable limits with no loss of information in the resultant en face projections; as a rule of thumb the axial resolution loss should not exceed the thickness of the thinnest useful en face projection -the thickness of ~20 µm is often used to generate images of the retinal capillary plexuses. In OCT imaging, specific terminology exists to describe the 3-D data sets. As in ultrasonography, a single image line obtained along the imaging beam, i.e. in depth, is known as an A-scan. It is represented as a 1-D matrix of pixels in the acquired data set. Traditionally, A-scans set the orientation of the Z-axis of the orthogonal coordinates system. A collection of A-scans acquired along a certain scanning line, for example along the fast scanning axis, forms a B-scan -a cross-sectional image. It is typically represented as a 2-D matrix of data pixels. The fasts scanning axis also typically sets the orientation of the X-axis. A set of B-scans acquired at subsequent locations in the object, typically along the slow scanning axis (Y-axis), forms a volumetric data set represented as a 3-D matrix. In OCTA, collections consisting of a few B-scans are typically acquired from the same location in the object to enable implementation of flow detection algorithms. Such groups of multiple B-scans are termed BM-scans. Cross-sectional images can be extracted from the 3-D data along any arbitrary slicing surface. If the surface is a plane perpendicular to the Z-axis (or depth-axis), then the resulting XY cross-section of single pixel thickness is typically called a C-scan. Stacks of C-scans selected from some arbitrary depth are used to generate en face projections. For clarity, the orientation of B-scans and C-scans in our data sets are depicted in Fig. 1 The data filtering pipeline is depicted in Fig. 2. The first two steps, Gaussian filtering of YZ sections and median filtering in Z-axis (a and b), are performed to remove the random noise, i.e. uncorrelated changes of the pixel values occurring among neighboring pixels. In practice, these initial steps remove the outliers of low and high intensity. The key filtering procedure is t reduce the no with the imag the speckle no character of the capillary v any given tim comparable w between the a the vessel. If  In the construction of the structuring filter, we use the knowledge that the orientation of the majority of retinal vessels is nearly perpendicular to the scanning beam of light. Thus, we apply the directional filtering in subsequent C-scans (XY planes indicated with pink rectangle in Fig. 1(b)). Using the knowledge that images of vessels form elongated shapes in the XYsections, we define an elliptical filter-kernel. The images of vessels can have all possible orientations in the XY planes, therefore we allow the kernel to rotate in a search for the most likely local orientation of the vessels in every point of the XY plane, and in all subsequent Cscans in the 3-D data set. The filtering is then done by calculation of the median within the kernels with the smallest value of an error function. Let C z (x, y) represent a C-scan at depth z in the volumetric OCTA data set. The orientation of our orthogonal reference system is indicated in Fig. 1. Here, we use discrete values of the x, y, z coordinates pointing to the voxels v(x, y, z).
We create a collection of elliptical kernels (2) n -an index numbering the ellipses in the collection, N -number of ellipses, a, c -the semi-major and semi-minor axes of the ellipse, β n -a rotation angle of the n th ellipse, relative to the x-axis, α -initial rotation angle.
The ellipses are discretized as illustrated in the inset of Fig. 2. All pixels fully encircled by the ellipse are given a value of 1. Pixels intersected by the ellipse are given a value of 1 if at least 50% of the pixel area falls inside the ellipse. All other pixels are given a value of 0. The number N of ellipses in the collection depends on the pixelization of the discretized kernel. The angles of rotation are selected in such a way that all the discretized ellipses contain a different collection of pixels, i.e. have different discretized shapes. In other words, the discretization of the rotation angle depends on the lengths (in pixels) of the ellipse axes a and c.
The collection of N ellipses is applied to every x, y point of the C-scan, in all subsequent C-scans in the 3-D data set. The ellipses are centered at the x, y points. To find the ellipse with the best match to the orientation of the local features in the OCTA images, we define a measure of similarity as an error function, , w h e r e : , , n n z n z n x y Err x y x y x y µ n , σ n -mean and standard deviation of the OCTA signal intensity within the n th discretized ellipse, σ z -standard deviation of the OCTA signal calculated from all pixels in the C-scan at depth z. The proposed metric uses local statistics to find the most likely match and uses global standard deviation value as a normalization term. When an ellipse lies along the vessel, the σ n (x, y) value is low and consequently the Err function has also a low value at the x, y point of the C-scan. For an ellipse that crosses the vessel, the value of σ n (x, y) is high and the value of the metric is closer to 1. The best matching ellipse, x y , at a given voxel v (x, y, z), is selected as the one with the smallest Err function value: where n z BEST (x, y) is the index of the best match in the collection of discretized ellipses {E 0 , E 1 , …, E N-1 }. The procedure of the best ellipse selection is illustrated in Fig. 2.
The final step of the elliptical median filtering is calculation of the median value of the OCTA signal within the best matching ellipse at all voxels, v (x, y, z), in subsequent C-scans, C z (x, y) , , where the curl brackets indicate a collection of voxels within the best matching ellipse at point x, y. After this step, the histogram of intensity in the three-dimensional data set is scaled to 0-255 range and the images are converted to 8-bit integer numerical-representation.

d) Generation of en face projection images
The next step of the data processing is generation of en face projection images. The surface approximating the shape of the outer plexiform layer was used to select depth ranges in the filtered 3-D OCTA data set. We generated angiograms at two depths as maximum intensity projections (MIP) of the OCTA data. The first visualizes the retinal vasculature and the inner capillary plexus. The second displays the outer capillary plexus. The outer plexiform layer segmentation lines and their depth positions used for generation of the en face projections are illustrated in Fig. 1(a, c).

e) Sharpening of the en face projections and contrast enhancement
In order to improve the definition of the vessels in the final en face projections, we enhanced high spatial frequencies (large image gradients) using a standard unsharp mask filter. First, the method uses image blurring by Gaussian smoothing. The blurred image is then subtracted from the original one to find the edges of the imaged features. The sharpened image is created by addition of the edge containing image, multiplied by empirically selected positive weight, to the original image. The thickness of the edges used for image sharpening is adjusted by tuning the width (standard deviation) of the Gaussian function. Then, a standard and commonly used contrast enhancement is implemented.

Practical implementation of the elliptical directional median filtering
We have implemented the elliptical directional median filtering using custom programming in C + + environment. The processing was performed automatically. In Table 3 we have listed values of filtering parameters used to generate the images presented in Figs. 4-6.  The ellipse parameters a and c were selected empirically. Several values were applied and the resultant images were presented to experts in the OCTA technique (MS, AS, IG) and an ophthalmologist (BS). The selection of the range of tested values was guided by the physical dimensions of the imaged vascular networks. The minor semi-axis a should be of an order of the smallest vascular lumen radius visualized in C-scans. The major semi-axis c can be estimated by taking into account the mean length of straight parts of the vessels and mean distance between the vascular junctions. Here, the selection was guided by visual inspection of the images. More advanced analysis of the images can be also performed to provide initial estimate of the tortuosity of the imaged vascular networks or to define the types, areas and distances between the vascular junctions. Note however that such analyses are typically preceded by the filtering of the images which is actually the goal of the image processing following these preparatory steps. A dilemma then arises how to perform the initial or recursive filtering.
Depending on the implemented scan protocols, the physical dimensions of the vascular networks translate into pixelization of their images. This pixelization determines the number N of rotated ellipses, i.e. the increment of the rotation angle. In practice, the guidance for the number N are the shapes of the discretized ellipses. On one hand, the collection should not contain identical kernel shapes since they will not influence the results but will increase the computational cost. On the other hand, the smallest difference of the kernel shapes should not be larger than one pixel. The latter criterion should be however tested in practice to determine if the single-pixel difference has a non-negligible impact on the filtering outcome. For the purpose of this paper, we have empirically selected the number N = 8 by visual comparison of the resultant images.

Filter-performance evaluation metric
To evaluate the performance of the elliptical directional filter depending on the kernel parameters, we propose to use a vascular continuity metric (VCont). We define VCont as a function which counts the number of vascular fragments. The vascular fragments are defined as clusters of connected pixels exceeding a selected intensity threshold value. With the increasing continuity of vessels, the number of vascular fragments will decrease. To measure whether the continuity of vessels has improved after application of the filtering procedure, the number of vascular fragments in the input image (A) is divided by the number of vascular fragments in the filtered image (B):

Demonstration of the elliptical directional filtering method in a normal eye
The effects of subsequent steps of the elliptical directional filtering method are summarized in Table 4. The summary shows the first of the two generated vascular layers, i.e. the angiograms of the retinal vasculature and the inner capillary plexus. The first row in the table illustrates the effects of filtering in an example single C-scan from the 3-D OCTA data set. Application of the first two preparatory steps, the YZ Gaussian filter and 1-D median filter, improves the visibility of vessels (image 1.b in Table 4). The "granular" noise is however still present within the vessels and in the image areas devoid of vasculature, e.g. in the foveal avascular zone forming a characteristic dark circle in the center of the angiograms. This can be attributed to minimal transverse filtering of the C-scans -the Gaussian filter affects only groups of 3 pixels in the Y direction (slow scanning direction). The effect of improved visibility of vessels can be attributed to the axial 1-D median filtering. It effectively combines information from stacks of neighboring C-scans but also degrades the axial resolution of the data set. Here, we have used a 9-pixel long window, corresponding to ~6 µm depth-range, which is comparable with the axial imaging resolution of our research-grade OCT imaging system. The key step of the processing pipeline is the elliptical directional filtering of the C-scans. It reduces the granular noise and smooths-out the images of vessels as demonstrated in angiogram 1.c. in Table 4.   Table 4. stribution Table 5

. Influence of the elliptical kernel parameter values on the OCTA en face projections. Parameters and evaluation-metrics values are given below the images.
Images obtained with the research grade OCT setup.

Reference image Symmetrical median filtering
No filtering c = a = 2 pixels c = a = 4 pixels The influence of various elliptical kernel parameters values on the output OCT angiograms is demonstrated in rows 2-4 of Table 5. For small values of c and a (e.g. 2.5 and 0.5 pixels), the results of the filtering are similar to the effects of application of a symmetrical median filter (c = a = 2 pixels). Increasing values of the c and a parameter improves the noise reduction. All values of the vascular continuity metric (VCont) are larger than 1, i.e. filtering reduces the fragmentation of vessels. In each row, VCont values increase with the growing parameter a. That is, increasing the width of the elliptical kernel with no change of its length, tends to increase the vascular continuity as compared to the unfiltered reference image. Here we tested values of a = 3.4 μm, 5.4 μm and 10 μm, comparable with the vascular radii of the retinal capillaries. On the other hand, in each column VCont values have a descending tendency, with the exception of the first column in which a local maximum exists for c = 26.8 μm and a = 3.4 μm. Intuitively, with increasing length of the kernel, filtering decreases the fragmentation of vessels. However, if parameter c becomes too large in combination with parameter a, then the images of vessels are excessively blurred and the visibility of the smallest vessels may decrease. In consequence, fragmentation of the vessels in the thresholded image increases causing a decrease of VCont values. The local maximum in the first column may suggest a good balance between the parameters c and a, resulting in a reasonable degree of filtering. In this case it also coincides with the parameters values experimentally chosen by the experts for the optimum filtering.  In unfiltered images, even well-defined vessels, e.g. insets (a) in Fig. 4 and 5, have noisy appearance. The elliptical directional filtering smooths-out the angiograms leaving the geometrical features of these vessels unaltered. Capillaries tend to appear fragmented in the unfiltered images. The fragmentation may appear as small, random breaks in otherwise well-defined vascular paths, e.g. insets (b) in Fig. 4 and Fig. 5, or they may appear as disconnected parts of the vasculature, e.g. insets (c) in Fig. 4 and 5. In the first case, our filtering method fills-in the gaps leaving all other characteristic features unchanged. In the second case the filtering finds the statistically most likely local orientations of vessels and generates the most probable shape correction, improving the continuity of vascular paths. The last demonstrated example are noisy regions with no apparent vascular networks in unfiltered angiograms, e.g. insets (d) in Fig. 4 and 5. Also here, the filtering results in generation of most probable fragments of vascular paths by statistical analysis of the OCTA intensity-distribution trends in the local neighborhood. Finally, it is worth noticing that in the foveal avascular zone (dark circular area in the center of the angiograms) no images of vessels are generated, which demonstrates that the elliptical directional filtering correctly differentiates between noisy vascular networks and noisy avascular areas in the OCTA data obtained in the imaging of normal eye.  Table 2). Case 3image size 5 x5 mm (scan protocol 2 in Table 2).

Imaging of diabetic retinopathy cases
The influence of the elliptical directional filtering on the OCTA imaging of retinal pathology is illustrated in Fig. 6. We show three diabetic retinopathy cases with increasing severity of vascular changes, imaged with Revo NX.
Case 1 is a 68 years old Caucasian female with non-proliferative diabetic retinopathy, early stage of the disease. We imaged the right eye, with best corrected visual acuity (BCVA) 0.9 (decimal scale). In this case filtering provides improved depiction of blood vessels and reveals the irregularity of the foveal avascular zone (FAZ). Visualizing the retinal capillary network around FAZ can be a useful tool in the very early diagnosis of diabetic retinopathy. For precise assessment of the vascular density and the FAZ area, the highest image quality with a good contrast is essential.
Case 2 is a 65 years old Caucasian male with non-proliferative diabetic retinopathy with enlarged foveal avascular zone. Images were taken in the right eye, BCVA 0.6. Reduced noise and improved continuity of vessels enables better assessment of the foveal avascular zone enlargement (pink arrow), a sign of increasing ischemia. Also, filtering helps in evaluation of the retinal capillaries' density -here decreased, as indicated with green arrows. The image quality became sufficient to reveal the presence of a few microaneurysms (blue arrows) in the nasal macula.
Case 3 is a 33 years old Caucasian male with proliferative diabetic retinopathy, with extensive capillary non-perfusion areas. Treatment with panretinal laser photocoagulation (PRP) was applied to the peripheral retina due to abnormal blood vessels. Images of the left eye with BCVA of 0.1 are shown. Here, the suppressed noise improves visualization of the non-perfusion regions (green arrows). The early detection of non-perfused areas in diabetic retinopathy is crucial because they are hallmark features of ischemic retinal disease. Improved continuity of the visualized vascular paths emphasizes the perfused but tortuous capillaries. This can be utilised for staging of the diabetic retinopathy and guiding the therapy.

Summary and conclusions
OCT angiography provides a set of challenges in the eye imaging. One of them is noise affecting the visualization of vessels, especially of small calibers -the capillaries. These random OCTA signal variations have an origin in the imaging technique but also in statistical character of motion of the blood cells through the narrow lumens of capillary vessels. We have provided an overview how the noise affects visualization of vessels in OCTA imaging of the eye fundus. For the demonstration we have selected the two-frame imaging regime. It minimizes the size of data sets, their acquisition times and enables maximization of the imaged areas, which are advantages desirable in clinical imaging. However, it is also one of the most noise-wise unfavorable OCTA scanning regimes and provides a challenge in image processing, which should reduce the noise without altering features relevant for the diagnostics. With these objectives in mind, we have developed a method of directional median filtering with elliptical structuring elements.
The median (or similar statistical operations, e.g. mean) typically provides the most effective filtering if a large number of points of the filtered image feature can be included to provide a good statistical sample of the signal. Intuitively, the smaller is the sample the less effective is the filtering. If the imaged features are elongated, then application of symmetrical kernels (e.g. square or circular) forces kernel sizes proportional to the smallest dimension of that feature. Otherwise, the information about vessels may be lost in the noise. Application of elongated kernels enables to include a larger number of meaningful data pixels and to avoid inclusion of meaningless noise in the filtering procedure. The advantage of using filters with kernel shapes matching the expected image features over symmetrical kernels is therefore in the inclusion of larger number of most likely useful data pixels.
We have demonstrated that our method reduces the noise and fragmentation of OCTA images by finding the most likely local orientations of vascular paths. In the areas where vasculature is present, the OCTA signal intensity typically shows local directionality which persists in the closest neighborhood. The rotation of ellipses enables finding those local trends in spatial distribution of the OCTA signal. In the straight sections of vessels, the elliptical kernels typically include only signal from the inside of the vessel in the median filtering. The influence of noise from the outside of the vessel is minimized. In the vascular junctions, the ellipses are aligned along the longest straight paths across the vessel branching. Since a part of the kernel may include the signal from the outside of the vessels, the filtering of the vascular junctions may be affected by the surrounding noise more than are straight parts of the vessels. If no trends exist, i.e. in the avascular regions of the image, the orientation of the ellipses changes randomly from voxel to voxel and returns a median of local noise distribution. As a result, the signal variations within the vessels are smoothed-out, vessels appear de-fragmented and the noise is tuned down.
To evaluate the influence of the elliptical kernel parameters on the resultant OCTA images, we have introduced a vascular continuity metric (VCont). With the use of this metric we have demonstrated that long and narrow kernels should be considered for median filtering of the OCT angiograms. However, care needs to be taken to avoid excessive lengths of the kernels which may result in "over-filtering" of the images and loss of the visibility of the small vessels. There can be of course further discussion whether the proposed metric is "good". For example, if its values reflect the quality of the images in a way which would be consistent with the judgement of an experienced human observer, or if it would enable finding the optimum filtering parameters automatically. The VCont metric measures only one possible feature of the filtered images -the continuity of vessels. While it correctly informs about increasing connectivity of vessels, it does not recognize when "overfiltering" occurs which may, for example, excessively blur the images of vessels, introduce erroneous connection between vessels or cause loss of visibility of small vessels. For this reason, VCont alone is not sufficient for the above-mentioned purposes. Additional metrics should be developed for application in combination with VCont to recognize the effects of over-filtering in the OCTA images. In such multi-parametric space extremum may exist pointing at the optimal filtering parameters. However, this topic is out of the scope of our paper and may be a subject of further study.
In the elliptical median filtering method, we have implemented only global contrast enhancement and aimed at preservation of the signal intensities gradation among the vessels. Vessels with strong OCTA signal in the unfiltered data have high pixel values in the filtered images. Vessels with weak OCTA signal prior to filtering have low pixel values after the filtering. Important information about changes in the blood flow may be included in the OCTA signal intensity. Therefore, we want to avoid excessive saturation of vessels in the images which may lead to loss of potentially useful information. Similarly, adjusting the contrast locally may change the relative differences in OCTA signal intensity among the imaged vessels and lead to the disruption of useful information. Moreover, image areas with weak signals are more likely suspects of erroneous information and therefore should be easy to identify and treated with greater caution in the interpretation of the results.
Erroneous information which may be introduced to OCTA images by any filtering procedure, including our elliptical median filtering method, is generation of false positive and/or false negative vessels, i.e. generation of images of vessels in the areas where no vessels are present in the imaged eye, or incorrect erasing of images of vessels. Such flawed identification of vessels is a potential limitation of our method which requires further investigation.
A typical issue in data filtering is degradation of the image resolution. In our method, the worst resolution loss is along the depth direction due to median filtering with a window length (6 µm) comparable to axial imaging resolution (5 µm). This however does not affect the visualization of retinal vascular networks in the en face projections. The resolution loss in the transverse direction (XY planes) depends on the selected parameters of the elliptical kernels. In the presented results, the minor axis of the ellipses was one pixel. The major axis was 8-pixels long but ellipses are aligned along the vessels. The loss of transverse image resolution can be minimized by appropriate selection of the ellipse parameters.
Numerically, elliptical filtering is straightforward to implement, although application of the elliptical filter should be performed in well motion-corrected data sets to provide the best results. I.e. images of the vascular networks should not contain discontinuities, for example caused by misalignment of the B-scans. If such eye motion induced artefacts exist in the data set, they will not be filtered out and will still appear in the filtered images as characteristic lines crossing the entire extent of the en face projections. The algorithm is easy to parallelize on central processing unit (CPU) and is well suited for execution on graphics processing unit (GPU). We performed our data processing on CPU with a code written in C++ and parallelized it with OpenMP programming interface. The implementation takes less than 10s to filter a 300x320x864 voxel data-matrix on a standard personal computer (PC) with i5-6600K CPU.
We believe that the simplicity of the method, easy implementation combined with decent results obtained in noisy data sets, makes the elliptical directional median filtering a useful tool in information-preserving quality improvement of the OCT angiograms.