Attenuation of stripe artifacts in optical coherence tomography images through wavelet-FFT filtering

The use of polarization-maintaining (PM) fibers for polarization-sensitive optical coherence tomography (PS-OCT) can result in numerous image artifacts which degrade the reliability of birefringence measurements. Similar artifacts can also arise in conventional OCT, due to stray reflections from optical surfaces, a problem which is increasing in tandem with the steady rise in source coherence lengths. Here, a recently presented wavelet-FFT filter[Opt. Express 17(10), 8567 (2009).19434191] is combined with surface flattening displacement fields in order to suppress ghost artifacts following either a duplicate or inverse profile to that of the sample surface. In addition, horizontal coherence stripes originating from Fresnel reflections of optical components are suppressed in order to facilitate accurate surface detection. The result is an improved visualization of the phase-retardance profile within tissue, which may improve the reliability of curve-fitting methods for localized birefringence estimation. While the results are presented with a focus towards PS-OCT, the filtering method can also be applied to the removal of stray reflection artifacts in conventional OCT images.


Introduction
Optical coherence tomography (OCT) is a non-invasive optical modality, capable of producing high-resolution cross-sectional images of tissue structure to a depth of a few millimeters [2]. One downside of regular, intensity-based OCT is a lack of tissue specific contrast, which often makes it challenging to discern between regions of tissue which exhibit similar refractive indices. As a result, numerous functional extensions to OCT have been developed, each aiming to exploit unique properties of the light and tissue in order to extract additional contrast. Such extensions include: angiographic OCT [3][4][5][6], spectroscopic OCT [7], OCT elastography [8] and polarization-sensitive OCT (PS-OCT) [9].
PS-OCT utilizes measurements of the polarimetric information carried by transverse light waves to measure properties such as the phase retardation (δ ) [10,11] and fast birefringent axis orientation (ϑ ) [12] within a sample. It has found extensive use within biological fields, notably within dermatology [13], opthalmology [14], dentistry [15] and for the enhanced study of tissues such as bone, tendon, ligament and cartilage [16]. Recently the use of singlemode (SM) fiber to construct PSOCT systems has increased in popularity [17][18][19], however PM fiber is highly stable in its birefringence properties and allows the relatively simple "single input state" method to be applied to fiber-based systems [20]. Indeed, many recently developed PS-OCT systems make use of polarization-maintaining (PM) fibers [20][21][22] which have less systematic errors and are often more clinically applicable when compared to systems utilizing bulk optics [23]. However such systems are not without disadvantage; the polarization-maintaining property of the fiber is achieved by inducing an intentional linear birefringence within the fiber, such that the two polarization modes (Horizontal and Vertical) propagate at distinct, and significantly different phase velocities. Effectively, beam propagation is slowed along one axis of the elliptic fiber core with respect to the other [21]. Slight imperfections in the splice angles and alignment of fiber connectors can result in undesirable cross-coupling between the two axis [21,22], which when detected and processed will result in vertically offset "ghost" copies of the OCT signal. While such "ghost" images are typically much lower in intensity that the original image, their bright surface profiles have sufficient intensity to noticeably degrade the underlying image. A hardware based method of removing such artifacts involves displacing them out of the OCT field-of-view through usage of long PM fiber segments [22]. Here, we propose a post-processing based algorithm based around wavelet-FFT filtering [1] which aims to attenuate the surface ghost reflections such that artifact-free images of the phase-retardance profile can be obtained, without the requirement for any hardware based compensation or modification. Fig. 1. illustrates the primary artifact types which this algorithm aims to address. All three artifacts are common occurrences within PM-fiber based PS-OCT images, with horizontal coherence noise stripes also affecting conventional OCT: 1. Horizontal coherence noise stripes -Uniform horizontal lines spanning the entire width of the image. Derived from Fresnel reflections from optical components within the system.

2.
Surface ghost reflections -Cross-coupling between the orthogonal polarization modes within the PM-fibers causes vertically offset copies of the OCT image to be present both above and below the sample. While these "ghost" images generally have a much lower signal intensity, their surface profiles are sufficiently bright such that the underlying image is significantly degraded. Thus the focus here is on the removal of the bright surface profile, rather than the entire "ghost" image.
3. Inverted surface ghost reflections -Complex-conjugate derived "mirror-artifacts" which arise from ghost artifacts which are positioned above the zero optical-pathdifference (OPD) point of the system. These become more problematic the closer the sample surface is to the location of the zero OPD.

PS-OCT system
A PM fiber based PS-OCT system was used to collect images of tissue for processing. The system follows the scheme originally reported by Al-Quasi et al [22] and has been detailed fully in a previous publication [24]. Figure 2 shows a schematic of the system. The light source used was a commercially available swept-source laser (HSL-2000-10-MDL, Santec Japan) which has a center wavelength of 1315nm, a sweep range of 157nm and a full width at half maximum (FWHM) of 128nm. The measured axial resolution (in-air) of the system was approximately 10μm and the laser sweeps at a rate of 10 kHz. Circularly polarized light is directed to the sample via a Panda PM-fiber Mach-Zehnder interferometer. Backscattered light from the sample is directed towards two balanced detection channels (1817-FC, New Focus, US) which measure the interferometric signal corresponding to the horizontal and vertical orthogonal polarization states of the light. In linearly birefringent materials, light with a linear polarization which is aligned parallel to the optic axis of the material experiences a different index of refraction (and thus propagation speed) compared to light with a linear polarization aligned perpendicular to the optical axis. Thus birefringence within the sample causes a phase-retardance between the complex amplitudes of each orthogonal polarization channel, which is detected by considering the evolution of the ratio between channels as a function of tissue depth. A MATLAB (R2015b -MathWorks) general user interface was designed to facilitate both data acquisition and live processing of the images with the algorithm discussed herein.

Wavelet-FFT filtering
Horizontal or vertical line or "stripe" artifacts are a common occurrence within many imaging modalities, including OCT. As a uniform stripe has a very high spatial frequency bandwidth along its width but a minimal frequency bandwidth along its length direction, a common method of removing such artifacts involves the use of the frequency domain. Thus, a basic implementation of a de-striping algorithm may simply calculate the 2D FFT ( , ) x y F u u of an image ( , ) f x y , then attenuate frequencies along the axis which is perpendicular to the line while maintaining the coefficients in the neighborhood of the origin. Either: The filtered image can then be recovered using the inverse FFT operation. One limitation of this approach is that image detail coefficients are prone to being deleted alongside the detrimental stripes, causing a loss of underlying image information. An improved method of filtering stripes from images was recently proposed by Münch et al [1]. They demonstrated that by first wavelet filtering the image, such that the approximation, diagonal detail, vertical detail and horizontal detail components of an image are reversibly condensed into separate bands, the aforementioned FFT filtering can be performed only on the relevant artifact corrupted band (horizontal in this case), thereby improving preservation of the underlying image information when compared to the purely FFT based method [1]. Figure 3 shows how this wavelet decomposition is performed. Fig. 3. A single-level wavelet decomposition of a PS-OCT structural image. Low and High refer to a Daubechies decomposition low-pass filter and highpass filter respectively, the image data is convolved with these filters along either the columns or rows. 1↓2 refers to a downsampling of the rows while 2↓1 refers to a down-sampling of the columns.

The image ( , )
f x y is fed through a bank of low and high pass filters which operate along either the rows or columns of the image. Each time the image is filtered, it is essentially down-sampled by a factor of two by the filter. Low pass filtering both the rows and columns captures the low frequency approximation information of the image. Low pass filtering the rows but high pass filtering the columns captures horizontal detail information, due to horizontal lines having a high frequency along the columns and a low frequency along the rows. Similarly, vertical detail information is extracted by high pass filtering the rows and low pass filtering the columns. Lastly diagonal information is extracted by high pass filtering both the rows and columns. This process can be repeated on the approximation band multiple times in order to extract lines which were wider than 1-pixel in the original image (Lines become a factor of 2 "thinner" between each successive wavelet decomposition due to the down-sampling). The aforementioned FFT filter can then be applied to only the relevant artifact corrupted band (Horizontal detail band here) and the filtered image repackaged using inverse wavelet transformations. This combined "wavelet-FFT" filter itself contains numerous variables: 1. The decomposition level ( l ) is the number of successive approximation bands which are wavelet filtered, a higher level will filter wider lines from the image.
2. Sigma ( σ ) defines the width of FFT coefficients which are damped within each horizontal detail band. A larger value of sigma will filter increasingly "imperfect" lines from the image, for example lines which are not uniform along their length. and a Daubechies 20 wavelet being selected as they offered a good balance between performance and processing speed. All horizontal detail bands were filtered in order to remove horizontal lines within the images. These parameters provided good results across a range of PS-OCT data sets, although improved results could potentially be obtained through optimization of the filter parameters for specific image cases. In addition, as the artifacts affecting the PS-OCT horizontal and vertical amplitude images are additive, they always appear brighter than the underlying true signal. Thus, a simple method of preventing the lines from "blurring" across previously unaffected signal [25] is to calculate the minimum value between the pre-and post-filtering images. This additional step was applied following each application of the wavelet-FFT filter. The benefits of this minimum filtering step are illustrated in Fig. 4 below. Fig. 4. A) A small section of a PS-OCT reflectivity image which is contaminated with horizontal stripe artifacts. B) The image following wavelet-FFT filtering, the filtered lines have "blurred" out into the surrounding signal. C) The result of taking the minimum value between A and B. D) Absolute difference image between A and B. Bright regions show where the image has changed. There is a visible "blurring" around each of the line artifacts which indicates that signal previously not corrupted by artifacts has been affected. E) Absolute difference image between A and C. Following the application of the minimum operation, the change in signal intensity is restricted to the location of the horizontal artifacts only.

Artifact removal algorithm
In order to efficiently filter the three artifact types shown in Fig. 1 without significantly degrading the underlying image, a multi-step image processing pipeline was applied to PS-OCT images derived from both the horizontal and vertically polarized channels. This processing pipeline is illustrated in Fig. 5 and described below. Firstly, the combined wavelet-FFT filter described in Sec. 2.2 is applied to both the horizontal and vertical polarized images ( Fig. 5A/D) in order to heavily suppress any horizontal coherence noise stripes which are present (Fig. 5B/E). The assumption at this point is that the sample surface within the image is not flat, otherwise it would also be removed by the filter. The next step involves detecting the sample surface, which is easier in the absence of the horizontal coherence noise stripes, however care must be taken to avoid erroneously detecting ghosted signal as the sample surface. The full sample surface detection process has been detailed previously [26], briefly this involves detecting the maximum signal intensity along each A-scan to generate a rough surface profile. This profile is then refined using a connectivity-enforcing path-finding algorithm. Once the surface profile is known, a displacement field ( D ) is defined as the pixel-mapping transformation which would flatten the sample surface. The inverse of this transformation ( 1 D − ) is first applied, flattening any inverted surface ghost reflections within the image, such artifacts are then attenuated through a second application of the wavelet-FFT filter. The now inverted surface is then flattened through two applications of the displacement field ( 2D ). While a third application of the wavelet-FFT filter at this point would indeed attenuate regular ghost artifacts within the image, it would also attenuate the true sample surface itself. To avoid this, the sample surface ± 5-pixels were stored prior to the wavelet-FFT filter being applied for a third time; the stored sample surface layer was then restored to its original location. The resulting horizontal (Fig. 5C) and vertical (Fig. 5F) images were now filtered of artifacts and could be processed into sample reflectivity (Fig. 5G) and phase-retardance (Fig.  5I) images.

Filter performance optimization
To optimize the performance of the filter, an OCT image acquired from a commercial system (VivoSight -Michelson Diagnostics Ltd) which did not present any noticeable image artifacts had artificial artifacts added to it. For this, horizontal coherence noise stripes were modelled as horizontal lines between 1 and 3 pixels in diameter with a random opacity of between 30 and 70%. Surface ghost reflections were simulated by offsetting copies of the B-scan by increments of 30-pixels in both vertical directions and setting the opacity randomly to between 10 and 90%. Inverse surface ghost reflections were produced in a similar manner, using an inverted B-scan. The opacity values were selected empirically based on the brightness of typical artifacts observed in real PS-OCT images. This artifact corrupted B-scan was then processed using the proposed filtering methodology before being compared with the original "clea a quality metr the underlyin 10 lo PSNR = between the p the tissue, no the original B trends toward the "clean" or a range of filt Decomposition level (l) The best p Daubechies 4 the PSNR of matches the o the artifact co for processed imag ise above the B-scan). In idea ds infinity, thu riginal. Table 1 er parameters. on of the ssed OCT ssing by a on of the value of parameters such as the decomposition level will vary depending on the spatial characteristics of the artifacts within the image (e.g artifact diameter).

Filter performance evaluation
To evaluate the performance of the algorithm using real PS-OCT data, a 6x6x2mm volume scan (Consisting of 1200x600x512 voxels) was collected from an ex-vivo sample of chicken breast tissue. Each slice of the volume was processed independently to generate both reflectivity and phase-retardance images following the methodology described in Sec 2.3. Figure 7 shows the results of this processing, to compare within the tissue only, the unfiltered phase-retardance images have been masked using the same mask as that used for the filtered phase-retardance images, this mask was generated as:  Within the unfiltered sample reflectivity PS-OCT image (Fig. 7A), bright horizontal coherence lines and ghosting artifacts are visible both within the sample and above the sample surface. Such lines increase the difficulty of robust sample surface detection and reduce structural clarity within the image. In addition, the artifacts are clearly visible within the unfiltered phase-retardance image (Fig. 7C), typically manifesting themselves as blue lines with a phase-retardance close to 0. The reason the artifacts appear blue is due to them having a comparatively higher brightness within the horizontal channel than the vertical (Fig. 5A/D), thus: Comparatively, both the filtered sample reflectivity image (Fig. 7B) and filtered phase retardance image (Fig. 7D) have visibly fewer artifacts than the unfiltered images, resulting in an improved visualization of the underlying tissue structure. The filtered images are not entirely artifact free, as small faint patches of ghost artifact can be seen towards the right side in addition to some faint banding at the upper left. From an en-face perspective across the entire volume, the artifacts are visible as wavy patches of bright reflectivity (Fig. 7E) or low phase retardance (Fig. 7F), due to them varying as a function of depth with respect to the sample surface. Following filtering, these wavy artifacts have been almost entirely attenuated (Fig. 7G/H) and are only visible as a faint blur in the reflectivity image.
Remaining artifacts which have persisted through the filtering process could potentially be addressed in two ways. Firstly, any minor imperfections in the sample surface detection will result in surface ghost reflections not aligning in a perfectly flat manner following the surface transformation process, resulting in non-uniform offset parameters in the frequency-domain. Thus, if further improvements can be made to accurately detect the sample surface, then the efficiency of filtering will likely increase. Secondly, further optimization of the wavelet-FFT filter parameters discussed in Sec. 2.2 could yield improved results, however the exact choice of parameters will likely depend on how far the various stripes affecting the images differ from ideal, 1-pixel wide horizontal lines. For example, increasingly wide line artifacts can be more effectively removed by increasing the decomposition level ( l ) of the filter at the expense of some degree of image information. A more efficient approach in terms of preserving the underlying image information, particularly for blurred lines such as those observed here, may be to specifically FFT filter the decomposition levels which correspond to the range of stripe widths within the image [1]. The entire filtering process, carries a processing overhead of ~1 second per 1200x512 pixel frame within MATLAB (Using a 3.5GHz i7 3770K CPU), thus this algorithm is unsuitable for high-speed real-time processing without further optimization; it is better suited for bulk offline processing of PS-OCT data sets following acquisition. In addition, although developed here specifically for filtering PMfiber based PS-OCT images, the algorithm could be useful for cleaning regular OCT images that are contaminated by stray reflections from straight and/or curved surfaces.

Conclusions
In summary, the proposed algorithm provides a means of filtering parasitic optical artifacts from OCT and PS-OCT data sets. The algorithm functions through a 3-step "stripe-removal" process whereby line or stripe artifacts are attenuated within an image using a combined wavelet-FFT filter. The image is filtered in its regular, inverted surface and surface flattened states to remove the different artifact types. Filtered images contained visibly attenuated coherence noise stripes together with a reduction in the intensity of both regular and inverted surface ghost reflections, with little apparent degradation to the underlying structural image. Importantly, the algorithm functions without any hardware modification to the OCT system, and can be easily applied retroactively to previously acquired scans. In general, this method enables a clearer visualization of both the structural sample reflectivity image and the corresponding phase-retardance plot. Curve-fitting methods used to estimate local birefringence within the tissue are likely to be more robust and require lower amounts of smoothing in the absence of the artifacts.
Future work will focus on the optimization of the filter parameters for specific data sets, allowing for more complete filtering of any present artifacts.