1. 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–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 ($\delta$) [10,11] and fast birefringent axis orientation ($\vartheta$) [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 single-mode (SM) fiber to construct PSOCT systems has increased in popularity [17–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–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:

 

Fig. 1 A) A structural PS-OCT scan of the palm showing the three common image artifacts. B) Artifacts of each type have been overlaid with a specific color: Green shows horizontal coherence noise stripes, red shows surface ghost reflections and blue shows inverted surface ghost reflections.

Download Full Size | PDF

2. Methods

2.1 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.

 

Fig. 2 A) Schematic diagram of the PM-fiber based PSOCT system. SS: Swept-source light source, PC: polarization controller, IL-LP: in-line linear polarizer, PMC1/2: polarization-maintaining couplers, QWP: quarter waveplate, PBS: polarization beamsplitter, H/V: Balanced photo-detectors which detect the horizontal and vertically polarized light respectively. B) MATLAB general user interface for acquisition and live processing of the data.

Download Full Size | PDF

2.2 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 $F(u_x,u_y)$ 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: $\text{{}F \text{|} u_x=0,u_y\ne 0\}=0$ for horizontal lines or $\text{{}F \text{|} u_x\ne 0,u_y=0\}=0$ for vertical lines. 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 down-sampling of the rows while 2↓1 refers to a down-sampling of the columns.

Download Full Size | PDF

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:

For this work, the wavelet-FFT filter parameters were optimized (See Sec. 3.1) with: $l=4$, $\sigma =10$ 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.

Download Full Size | PDF

2.3 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 ($D^{-1}$) 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 $\pm$ 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.

 

Fig. 5 The processing pipeline for PS-OCT artifact reduction. A) Unprocessed PS-OCT scan of a fingertip showing horizontally polarized signal ($A_H$). B) $A_H$ following the first application of the wavelet-FFT filter, which aims to attenuate horizontal coherence noise stripes. C) $A_H$ following the second and third applications of the wavelet-FFT filter, which aim to attenuate both regular and inverted surface ghost reflections. Signal located above the detected skin surface has been set to zero. D-F) Same as A-C but showing the vertically polarized signal ($A_V$). G) Processed sample reflectivity image calculated as: $R=A_H^2+A_V^2$. H) Unprocessed phase-retardance map calculated as: $\delta ={\tan }^{-1}(A_V/A_H)$. I) Processed phase-retardance map.

Download Full Size | PDF

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.

3. Results and discussion

3.1 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 “clean” B-scan through calculation of the peak signal-to-noise ratio (PSNR). PSNR is a quality metric which is commonly used to evaluate image processing algorithms, providing the underlying image does not change between comparisons [27]. It is calculated as: $PSNR=10{\log }_{10}(\frac{{255}^2}{MSE})$ for the 8-bit OCT images, where MSE is the mean square error between the processed image and the original B-scan. To enable comparison primarily within the tissue, noise above the skin surface was set to zero for all PSNR calculations (Including the original B-scan). In ideal cases, where the mean square error trends towards zero, PSNR trends towards infinity, thus higher values generally mean the image more closely matches the “clean” original. Table 1 shows the result of calculating PSNR for images processed with a range of filter parameters.

Table 1. Optimization of the filter parameters (, and wavelet type)

View Table

The best performance was observed following filtering to 4 decomposition levels, using a Daubechies 40 wavelet with $\sigma =10$, this corresponded to a PSNR of 36.64. For reference, the PSNR of the unprocessed image was 27.19 thus the processed image more closely matches the original image. Figure 6 shows the visual difference between the original image, the artifact corrupted image, and the image processed using above parameters.

 

Fig. 6 A) B-scan image of skin absent of any visible artifacts. B) The same image as A but artificial artifacts have been added to closely resemble a corrupted PS-OCT B-scan. C) The result of processing B with filter parameters: $\sigma =10$, $l=4$ and a Daubechies 40 wavelet.

Download Full Size | PDF

One aspect not discussed above is processing time, which varies as a function of the parameters. With this in mind, a Daubechies 20 wavelet was used for all other processed OCT images in this manuscript as it offers comparable PSNR (36.25) but speeds up processing by a factor of ~2x when compared to a Daubechies 40 wavelet. In addition, a limitation of the above analysis is the usage of only a single image. It is likely that the optimal 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).

3.2 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: $R_{Filtered}\left(x,y\right)>0.05$. The value of $0.05$ was measured as the noise-floor by calculating the 95th percentile of pixel brightness in the bottom 1200x600x50 pixels of the data set.

 

Fig. 7 Results of processing PS-OCT images acquired from ex-vivo chicken breast tissue. A) Unfiltered sample reflectivity image. B) Filtered sample reflectivity image. Reflectivity calculated as: $R=A_H^2+A_V^2$. C) Unfiltered phase-retardance image, for comparison this has been masked using the same mask as that used for D. D) Filtered phase-retardance image, masked by excluding areas of low signal from the image (Defined as a processed sample reflectivity of less than 0.05). Phase retardance calculated as: $\delta = {\tan }^{-1}(A_V/A_H)$. E) En-face view of the unfiltered reflectivity volume at a depth 55-pixels below the sample surface. F) En-face view of the unfiltered phase retardance volume at a depth 55-pixels below the sample surface. G) En-face view of the filtered reflectivity volume at a depth 55-pixels below the sample surface. H) En-face view of the filtered phase retardance volume at a depth 55-pixels below the sample surface. The pink vertical line on E/G shows the B-scan location used for the cross-sectional images.

Download Full Size | PDF

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: ${\tan }^{-1}(A_V/ A_H)\to 0$. 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 PM-fiber 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.

4. 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.