Open Access
Add to BibSonomyAbstract
We present a high speed polarization sensitive spectral domain optical coherence tomography system based on polarization maintaining fibers and two high speed CMOS line scan cameras capable of retinal imaging with up to 128 k A-lines/s. This high imaging speed strongly reduces motion artifacts and therefore averaging of several B-scans is possible, which strongly reduces speckle noise and improves image quality. We present several methods for averaging retardation and optic axis orientation, the best one providing a 5 fold noise reduction. Furthermore, a novel scheme of calculating images of degree of polarization uniformity is presented. We quantitatively compare the noise reduction depending on the number of averaged frames and discuss the limits of frame numbers that can usefully be averaged.
©2011 Optical Society of America
1. Introduction
Optical coherence tomography (OCT) has developed into a very useful technique for imaging of transparent and scattering tissue [1, 2], and can now be regarded as a standard tool in retinal diagnostics [3]. The main advantages of OCT are its high sensitivity, high acquisition speed and its high resolution down to a few microns. Spectral domain (SD) OCT techniques [4–6] have enabled high speed imaging in 3 dimensions with several 100k A-lines per second [7].
Polarization sensitive (PS) OCT is a functional extension of OCT [8, 9, 10]. It takes advantage of the additional polarization information carried by the reflected light, and can therefore add new image contrast compared to intensity based OCT. PS-OCT can reveal important information about biological tissue, such as quantitative distribution of birefringence, which is unavailable in conventional OCT.
A very interesting application field for PS-OCT is retinal imaging [11, 12]. One can distinguish between polarization preserving (e.g. photoreceptor layer), birefringent (e.g. retinal nerve fiber layer (RNFL), Henle’s fiber layer) and depolarizing layers (e. g. retinal pigment epithelium (RPE)). In addition, PS-OCT can provide qualitative and quantitative information on birefringent [13–15] and depolarizing tissues [16–21].
First PS-OCT systems were based on bulk optics [8–10], which is not as convenient in terms of easy alignment and handling compared to fiber based systems. Therefore different fiber based PS-OCT system have been proposed [22–24].
Recently time and spectral domain (SD) PS-OCT systems based on polarization maintaining fibers have been presented [25–27]. PM fiber based PS-OCT allows to maintain one of the main advantages of previous bulk-optics PS-OCT [17], i.e. the use of only a single input polarization state to simultaneously acquire reflectivity, retardation, optic axis orientation, and Stokes vector. Various methods have been proposed to correct the phase shift between the orthogonal polarization states traveling in the two orthogonal modes of the PM fiber (e.g. insertion of birefringent wedges, exact length matching of reference arm and object arm fiber). In our recent work, this phase shift is compensated by software in post processing.
One of the main limiting factors of PS-OCT image quality and quantitative data is speckle noise. Various methods of speckle noise reduction, based on hardware or software methods have been proposed for intensity based OCT [28–33]. A very successful method is averaging of several B-scans that are successively recorded at the almost the same position. This method is implemented in a commercially available OCT retinal scanner (Spectralis, Heidelberg Engineering) [34]. By employing a retinal tracker, the instrument is able to record several successive B-scans almost on the same retinal location. Typically 20 B-scans are averaged to provide a high-quality nearly speckle free image of the retina.
In this paper, we demonstrate a related method for generation of high quality speckle reduced PS-OCT images. Instead of using a hardware eye-tracker, we take advantage of new high-speed CMOS line scan cameras. We record images at a scanning speed of 70-128 kAlines/s, up to more than 6 times faster than with our previous systems for retinal imaging. This reduces motion artifacts considerably. Residual motion artifacts are corrected by a software algorithm. Up to 60 successively recorded and motion corrected B-scan data sets are averaged. We record images at nominally the same position (i.e. at similar settings of the x-y galvo scanner). Due to eye motion, the speckle pattern of each B-scan changes slightly and therefore after averaging one can obtain speckle reduced PS-OCT images of intensity, retardation, and optic axis orientation. Different averaging schemes are compared, the best ones providing a ~5 fold noise reduction in polarization sensitive images. Furthermore, a novel scheme of calculating images of degree of polarization uniformity (DOPU), is presented. Instead of spatial averaging of Stokes vector data, temporal averaging over successive images is used, which provides a considerably improved spatial resolution of DOPU images. Finally the limits of these methods are discussed.
2. Methods
The modified system is based on the instrument presented in [26], with the main difference of a modified spectrometer design including new CMOS cameras. We used a light source with a bandwidth of 50 nm and a center wavelength at 840 nm. Light emitted from the source enters a fiber based isolator to prevent damage of the source due to backreflection. A polarization control paddle is implemented to match the polarization state of light with the orientation of the fiber based polarizer. The polarizer provides vertically polarized light which enters the PM-fiber based interferometer, where it is split by a 50/50 coupler into a reference and a sample arm. In the reference arm, light exits the PM fiber via a collimator and passes a variable neutral density filter (for optimizing reference power), a quarter wave plate (QWP) oriented at 22.5°, a pair of glass prisms for variable dispersion compensation, and is finally reflected by the reference mirror. After double passage of the QWP, the orientation of the polarization plane is at 45° to the horizontal, providing equal reference power in both channels of the polarization sensitive detection unit. In the sample arm, after exiting the fiber through a collimator, light passes a QWP oriented at 45° which provides circularly polarized light onto the sample. With the help of an x-y galvanometer scanner and a telescope the beam is scanned over the retina. After recombination of light from the reference and the sample arm at the 2x2 fiber coupler, light is directed via a PM fiber to a fiber based polarizing beam splitter. The two orthogonally polarized beams are guided into two separate, identical spectrometers. Each spectrometer consists of a transmission grating (1200 lines/mm), a camera lens with a focal length of 150 mm, and a 4096 element line scan CMOS camera (Basler sprint spl 4096). To reduce the data transfer rate, only 2048 pixels of each camera were read out. Our system was operated at an A-scan rate of 70 k A-lines /sec and 128k A-lines/sec, covering a scan field of 15°x15°. One B-scan consists of 1000 A-scans (14 ms acquisition time @ 70kHz and 7.8ms @128kHz). For speckle reduced images, 60 B-scans are averaged (840 ms acquisition time @ 70 kHz and 468ms @128kHz), and 3D data sets consist of 120 B-scans (acquisition time 936ms). The sensitivity of the system was 96 dB at 70k A-lines/sec and 94 dB at 128k A-lines/sec. The depth resolution of the system is 4.5µm in tissue and the transversal resolution is 12µm.
3. Signal Processing
After standard SD-OCT data preprocessing (fixed pattern noise removal, zero padding, rescaling of the spectral data from wavelength to wave number space, dispersion compensation) an inverse FFT is performed. In a next step, the PM fiber length mismatch is compensated (a detailed description of this method can be found in [26]), and from the resulting structure term the amplitude and phase values of each channel are calculated. From these quantities the reflectivity, the retardation, and the optic axis orientation can be calculated as described in [35].
With the system used for this work an imaging speed of up to 128 kHz A-scans/s was achieved. This high acquisition speed enables imaging the retina with minimal motion artifacts, which allows to average several B-scans recorded at nominally the same location, which strongly improves the image quality by increasing the signal to noise ratio and reducing speckle noise. Therefore, up to 60 B-scans where acquired at the same position and averaged afterwards. To minimize remaining motion artifacts custom made correction software was applied [36]. The software is based on the calculation of a cross correlation of each individual B-scan to a fixed reference frame. This fixed reference frame can be arbitrarily chosen from the stack of B-scans. First, the amplitude square sum of channel 1 and 2 is calculated for each B-scan afterwards the resulting images are motion corrected. The resulting translation matrix for each B-scan is saved. With the help of these matrices all other parameters can be motion corrected (i.e.: amplitude of channel 1 and channel 2, retardation, optic axis orientation, and each Stokes vector element). After motion correction, the image stacks are averaged. While averaging of the intensity data is straightforward, there are different possibilities to average the polarization data which yield different image quality. First, different averaging techniques for retardation are discussed. The first method is based on calculating the retardation for each individual B-scan and averaging all calculated images afterwards (retardation averaging method 1)
were AV (x, z) denotes the amplitude of the vertical polarization channel, x the transversal position, z the depth position, AH (x, z) denotes the amplitude of the horizontal polarization channel, i is the frame number and n is the number of averaged B-scans. The second method averages all images of the vertical (AV) and horizontal (AH) channel separately. Afterwards one can calculate the retardation δ with the following formula (retardation averaging method 2):The third method is to apply a histogram based algorithm (similar to modal averaging). First, in each individual retardation B-scan the retardation at each pixel coordinate δi (x, z) (i…frame number) is taken and stored in a separate linear array Aδ (x, z) holding all retardation values of the corresponding position (x,z). Afterwards for each array Aδ a histogram of retardation values is calculated. The peaks of the histograms obtained for each array Aδ are regarded as the corresponding retardation value (retardation averaging method 3).With the setup used in this work the optic axis orientation θ can be calculated with the following formula
were ΔΦ is the phase difference of the vertical and horizontal channel. As can be seen from the above equation, the optic axis orientation is entirely encoded in the phase difference. The 2π ambiguity of the phase difference of the two channels makes simple averaging not meaningful. Due to noise, phase values near the limits -π or π can randomly jump from one of the extremes to the other and therefore an averaging would yield values near 0. To solve that problem one can either use a similar histogram based algorithm as described above for averaging retardation (optic axis averaging method 1), or one can use an extended averaging algorithm for phase values as described for Doppler imaging [37]. For this method a complex phasor (x,y) for each phase difference is calculated:Afterwards the average phase difference ΔΦAV(x,y) can be calculated with the following formula (optic axis averaging method 2)This method is similar to the Kasai estimator method used for Doppler OCT imaging [38], however applied to unit amplitude phasors.As an additional polarization parameter we calculated the degree of polarization uniformity (DOPU) [19], a parameter closely related to the degree of polarization (DOP) (which cannot be directly measured by a coherent detection technique like OCT). According to the original definition, DOPU is calculated from Stokes vector elements averaged within a tiny evaluation window which floats over the whole B-scan. This value can be used to segment depolarizing layers such as the retinal pigment epithelium. A detailed description of this method can be found in [19]. This method has the drawback of reduced spatial resolution of DOPU images because the images are a convolution of real structure with an evaluation window size. Here we chose a different approach to calculate DOPU maps with improved resolution. Since our new high speed PS-OCT setup allows very fast acquisition of B-scans at nominally the same location, we can replace the spatial evaluation window for calculating DOPU with a temporal window. Due to eye motion between successively acquired B-scans the speckle pattern at each xz position changes slightly and therefore a meaningful ensemble of speckles can be used for calculating DOPU. Similar as in the histogram based algorithms described above, arrays Si (x, z) of Stokes vectors are built. From these, sub-arrays of normalized Stokes vector elements Qi(x,z), Ui(x,z),Vi(x,z) are calculated. DOPU (x,z) is now calculated by:
Were the index m indicates that the respective Stokes vector element was averaged within the respective sub array (i.e. the averaging of Stokes vector elements is performed over pixels recorded successively at nominally the same spatial position). This enables to calculate DOPU images with improved spatial resolution. The actual resolution depends on the residual motion artifacts after software motion correction. It should be mentioned that there is a trade-off between actual resolution and number of independent speckles included in the averaging process (see discussion).An interesting parameter for glaucoma diagnosis is birefringence of the RNFL. The typical approach for calculation of the birefringence is described in [13–15]. Briefly, first one has to align each A-scan of the retardation image with respect to the detected retinal surface in a way that the RNFL surface appears as a horizontal line in the B-scan. After this alignment, the retardation data are averaged in transverse direction along the B-scan (floating average over typically 75 adjacent A-scans) to reduce speckle noise. Then a linear regression in depth direction is performed on the data points between surface and posterior border of the RNFL. The slope of the linear regression is proportional to the birefringence value at this location. Since our methods described above provide already reduced speckle noise by averaging the intensity and retardation data, no further averaging for calculation of the birefringence is necessary, and therefore an A-scan wise calculation of birefringence is possible, thus improving the spatial resolution of birefringence measurement.
4. In Vivo Imaging
PS-OCT imaging was performed in a healthy volunteer after informed consent was obtained. All measurements were approved by the ethics committee of the Medical University of Vienna and followed the tenets of the Declaration of Helsinki.
A promising application of PS-OCT is glaucoma diagnosis. In glaucoma, the RNFL is damaged, leading to reduced RNFL thickness and birefringence. It is hypothesized that a decrease of birefringence precedes a thickness reduction of the RNFL [39]. This effect can best be observed in areas of thick RNFL bundles around the optic nerve head. Therefore, circular scans around the optic nerve head are standard protocols of glaucoma diagnostics by OCT.
Figure 1 shows a circular scan around the optic nerve head (healthy subject) with a diameter of ~10 deg (~3 mm diameter). The images were acquired with 70k A-lines per second. Figure 1(a) shows the intensity image of a single circular scan, and Fig. 1(b) shows the intensity image averaged over 60 B-scans.
Fig. 1 Circumpapillary PS-OCT scans (1000 A-scans from healthy human retina in vivo. Intensity image on log scale (a) single scan; (b) averaged over 60 circular scans; Orientation of scan from left to right: Inferior, nasal, superior, temporal, inferior. (c) zoomed in region of Fig. 1(a) (rectangular window), (d) zoomed in region of Fig. 1(b) (rectangular window). Scan diameter: ~10 deg (corresponds to a circumference of ~9.4 mm, equal to horizontal image width; optical image depth: 1.8 mm).
The image shows increased thickness of the RNFL, the topmost bright reflecting layer, in the superior and inferior region, and a thin RNFL in the nasal and temporal region. Also several blood vessels can be observed. Figure 1(b) shows that averaging strongly improves the signal to noise ratio and also reduces speckle noise. Figure 1(c), 1(d) show a zoomed region of Fig. 1(a) and 1(b) (white rectangular window). The ghost images in the area in front of the retina are probably caused by imperfect optical elements that cause cross coupling of polarization states into the wrong mode of the PM fibers. A solution for removing these ghost images would be using longer PM fibers in the Michelson interferometer which would further separate the position of the ghost images (for a detailed discussion see [26]).
Figure 2(a) shows a single circular scan of retardation (same data set as in Fig. 1). To avoid erroneous data points of retardation and optic axis orientation caused by noise we first apply a thresholding procedure based on the intensity data to gate out areas with low signal intensity (threshold value: 3 times local intensity noise). Areas below the intensity threshold are displayed in grey. To reveal the true birefringent properties of the RNFL we had to compensate for corneal birefringence. We used the retardation and axis orientation measured locally at the retinal surface around the nerve head to correct for anterior segment birefringence by a software algorithm [40]. The expected increase of retardation with depth at the thickest RNFL bundles can be observed in the retardation image. To compare the performance of the different averaging techniques described in section 3 we applied the different averaging algorithms to the data set. Figures 2(b)–2(d) show images where retardation averaging methods 1, 2 and 3 have been used. One can clearly see that all three averaging techniques strongly improve the overall image quality, compared to the single circular scan. Figure 2(e), 2(f) show a zoomed region of Fig. 2(a) and 2(c), respectively (black rectangular window).
Fig. 2 Circumpapillary retardation PS-OCT scans. (a) retardation of single circular scan; (b) retardation averaging method 1; (c) retardation averaging method 2; (d) retardation averaging method 3; color bar for all scans 0°-90°.) (e) zoomed in region of Fig. 2(a) (rectangular window), (f) zoomed in region of Fig. 2(c) (rectangular window). Rectangular window in (b) (not drawn to scale) indicates position where statistics was performed (see text). Scan diameter: ~10 deg (corresponds to a circumference of ~9.4 mm, equal to horizontal image width; optical image depth: 1.8 mm).
Figure 3(a) shows a single circular scan of optic axis orientation (same data set as in Figs. 1 and 2). The optic axis orientation image shows two full color oscillations from left to right, in good agreement with the radial orientation of the nerve fiber bundles around the optic nerve head (corresponding to 2 x 180° axis orientation change). For the optic axis orientation data we also applied the different averaging techniques described in section 3 and compared the results. Figures 3(b), 3(c) show images where optic axis orientation averaging methods 1, 2 have been used. Figure 3(d) shows the same data as Fig. 3(c) however, with a cyclic color scale [41]. An advantage of such a cyclic color scale is that there are no abrupt color jumps at the -π/+π border; instead, a smooth change is observed, in agreement with the continuous transition along the circle circumference. Figures 3(e), 3(f) show a zoomed region of Fig. 3(a) and 3(d), respectively (black rectangular window).
Fig. 3 Circumpapillary optic axis orientation PS-OCT scans. (a) optic axis orientation of single circular scan; (b) optic axis orientation averaging method 1; (c) optic axis orientation averaging method 2; (d) optic axis orientation averaging method 2 with cyclic color bar (color bar on bottom right); color bar for all scans −90° to +90° (e) zoomed in region of Fig. 3(a) (rectangular window), (f) zoomed in region of Fig. 3(d) (rectangular window). Scan diameter: ~10 deg (corresponds to a circumference of ~9.4 mm, equal to horizontal image width; optical image depth: 1.8 mm).
We also quantitatively compared the signal improvement achieved by the averaging techniques. For this purpose we cut out a window with a size of 50(x-direction) x 100(z-direction) pixels at a position where the retardation and optic axis orientation should stay rather constant (see white rectangle in Fig. 2(a), window not drawn to scale) and calculated the mean value and the standard deviation of these values within the evaluation window. This procedure was applied for each image in Fig. 2 and 3.
The results for the retardation averaging methods are shown in Table 1 . As can be seen all averaging methods strongly decrease the standard deviation. Retardation averaging method 2 yields the best results and reduces the standard deviation by a factor of 5. Therefore for the remaining part of this paper only retardation averaging method 2 was used. Similar analysis was performed for the optic axis orientation averaging methods. The results are shown in Table 2 . It can be seen that optic axis orientation averaging method 2 yields the best results and decreases the standard deviation by a factor of ~3.5. Therefore for the remaining part of this paper only optic axis orientation averaging method 2 was used.We also quantitatively compared the signal improvement depending on the number of averaged frames. For this purpose a similar analysis as described above was performed. The results are shown in Table 3 and are discussed in detail in section 5. Figure 4(a) , 4(c), 4(e) shows images of intensity, retardation, and optic axis orientation averaged over 40 frames and Fig. 4(b), 4(d), 4(f) images averaged over 10 frames (same data set as in Figs. 2 and 3). As can be seen, the quality of the images averaged over 40 frames looks similar to the images averaged over 60 frames, whereas when only averaging over 10 frames the image quality is considerably decreased.
Table 1. Quantitative Comparison of Retardation Averaging Techniques
Table 2. Quantitative Comparison of Optic Axis Orientation Averaging Techniques
Table 3. Quantitative Comparison of Number of Averaged Frames
Fig. 4 Circumpapillary PS-OCT scans. (a), (c), (e):Intensity, retardation, optic axis orientation averaged over 40 circular scans; (b), (d), (f): Intensity, retardation, optic axis orientation averaged over 10 circular scans. Scan diameter: ~10 deg (corresponds to a circumference of ~9.4 mm, equal to horizontal image width; optical image depth: 1.8 mm).
As mentioned before, birefringence, i.e. retardation per unit depth, is an important quantity that might be useful for glaucoma diagnostics and that can be measured directly by PS-OCT. Figure 5 shows birefringence evaluation of the RNFL along the circumpapillary scan of Figs. 1, 2, and 3. The boundaries of the RNFL which were detected by an intensity based algorithm (as described in [15]) are indicated by red lines (Fig. 5(a)). Figure 5(b) shows a circumpapillary plot of birefringence without any lateral averaging.
Fig. 5 (a) Circumpapillary intensity scan (averaged) with segmented RNFL (red lines); (b) Circumpapillary plot of birefringence.
Another important region for retinal diagnostics is the foveal area. Figure 6 shows images of the fovea of a healthy volunteer. Figure 6(a) shows the intensity image of a single B-scan scan, and Fig. 6(b) shows the intensity averaged over 60 B-scans.
Fig. 6 Intensity image of the fovea of a healthy volunteer @ 70 kHz A-scan rate. (a) Single scan; (b) Intensity averaged over 60 B-scans. Image size: 15° (horizontal) x 0.75 mm (vertical, optical distance).
In the single scan image (6a) all the layers and features known from non polarization sensitive ultra high resolution OCT systems can be observed. Figure 6(b) shows that averaging strongly improves the signal to noise ratio in general and especially reduces speckle noise.
Figure 7 shows images of retardation and optic axis orientation obtained from the same data set as used in Fig. 6. Figures 7(a), 7(c) show the retardation and optic axis orientation of a single B-scan, respectively. Figure 7(b) (retardation) and Fig. 7(d) (optic axis orientation) show images averaged over 60 B-scans. The single scan images show, that most of the retinal layers in this area do not alter the polarization state. However, the RPE changes the polarization state of backscattered light in a random way (random variation from speckle to speckle), indicating depolarization [16]. Due to speckle noise reduction the averaged images of retardation and optic axis orientation look much smoother. However, due to the averaging process the depolarizing effect of the RPE is washed out.
Fig. 7 PS-OCT B-scan image of healthy human fovea in vivo. (a) retardation of single B-scan; (b) retardation averaging method 2; (c) optic axis orientation of single B-scan; (d) optic axis orientation averaging method 2; color bars similar to Figs. 2 and 3. Image size: 15° (horizontal) x 0.75 mm (vertical, optical distance).
From the PS-OCT data the DOPU image can be calculated [19]. The DOPU image quantitatively shows the extent of depolarization caused by retinal layers, and was shown to be a helpful tool for segmenting the RPE. Figure 8(a) shows a DOPU image of a healthy human retina obtained by the originally reported method from a single B-scan set of Stokes vectors with an evaluation window size of ~70(x) µm x 18µm(z). The depolarizing RPE with its low DOPU value (green color) is clearly observed. Figure 8(b) shows a DOPU image which was calculated from 60 successive B-scans as described in section 3. Since in this case no spatial averaging within one B-scan is necessary, the resolution of the resulting DOPU image is strongly improved. Figures 8(e), 8(f) show a zoomed region of Fig. 8(a) and 8(c), respectively (black rectangular window).
Fig. 8 PS-OCT B-scan of healthy human fovea in vivo. (a) degree of polarization uniformity DOPU (color bar: 0-1) calculated from single frame (b) DOPU image calculated from 60 B-scans. (c) zoomed in region of Fig. 8(a) (rectangular window), (d) zoomed in region of Fig. 8(b) (rectangular window). Image size: 15° (horizontal) x 0.75 mm (vertical, optical distance).
To further reduce remaining motion artifacts we increased the imaging speed to 128k A-scans per second. Figure 9 shows intensity images of the fovea of a healthy volunteer. Figures 9(a) and 9(b) show a single B-scan scan, and an average of 60 B-scans, respectively.
Fig. 9 Intensity image of the fovea of a healthy volunteer @ 128 kHz A-scan rate. (a) single scan ; (b) average of 60 B-scans. Image size: 15° (horizontal) x 0.75 mm (vertical, optical distance).
Although S/N is reduced at 128 kHz scan rate, still all the major layers are visible. Figure 9(b) shows that the averaged image is nearly without any motion artifacts. Compared to 6b the higher imaging speed seems to improve the visibility of fine structures like capillaries in the anterior retinal layers, probably because less residual motion artifacts are present.
Figure 10 shows polarization sensitive images retrieved from the data set of Fig. 9. Figures 10(a) and 10(c) show retardation and optic axis orientation of a single B-scan. The polarization scrambling caused by the RPE is still visible, although S/N is reduced. Figures 10(b) (retardation) and 10d (optic axis orientation) are averaged over 60 B-scans. Figure 10(e) shows a DOPU image calculated from 60 successive B-scans.
Fig. 10 PS-OCT B-scan image of healthy human fovea in vivo @ 128 kHz A-scan rate (a) Retardation of single B-scan; (b) retardation averaging method 2; (c) optic axis orientation of single B-scan; (d) optic axis orientation averaging method 2; (e) DOPU image calculated from 60 B-scans; color bars similar to Figs. 2–3. Image size: 15° (horizontal) x 0.75 mm (vertical, optical distance).
To demonstrate the 3D imaging capability of our instrument we imaged the optic nerve head region of a healthy volunteer. The 3D image stack consists of 120 B-scans (1000 A-scans per B-scan) and was recorded at an A-scan rate of 128 kHz, which results in 930ms for the whole 3D data set. A reconstructed fundus intensity image of the optic nerve head region can be seen in Fig. 11(a) . The fundus image is obtained via intensity summation along A-scans (z-projection) [42].
Fig. 11 Comparison of phase retardation introduced by the RNFL in the nerve head region of a healthy volunteer measured with SLP and PS-OCT (a) Fundus image reconstructed from the 3D data set; (b) Phase retardation map calculated from the PS-OCT data (color bar: deg); (c) Phase retardation map recorded with SLP (color bar: deg).
From this 3D data set we also calculated phase retardation maps showing retardation introduced by the RNFL. A detailed description of this method can be found in [43]. We further compared the results with scanning laser polarimetry (SLP) [44], a commercially available SLO based polarimetry technique (Carl Zeiss GDx VCC). Figures 11(b) and 11(c) show phase retardation maps obtained by PS-OCT and SLP in the same eye. In the SLP image (Fig. 11(c)) the expected nerve fiber layer retardation pattern can be observed (i.e. increased retardation superior and inferior to the nerve head, indicating thicker RNFL compared with decreased retardation temporally and nasally, indicating thinner RNFL). In the PS-OCT image (Fig. 11(b)) a very similar pattern can be observed and both images are in good agreement.
5. Discussion
We have presented an improved SD PS-OCT system capable of retinal imaging with up to 128 k A-lines/s. The improved imaging speed reduces motion artifacts and enables to record several ten B-scans at nominally the same location within a fraction of a second. Residual motion artifacts can be compensated by software postprocessing. The acquisition of several B-scans at the same position allows to generate averaged images with greatly improved image quality which is due to reduced noise. While noise reduction by averaging is already commercially available for intensity based retinal imaging (Spectralis OCT, Heidelberg Engineering), it has not yet been demonstrated for PS-OCT. In this paper, we demonstrated noise reduction by averaging for PS-OCT and, furthermore, compared different averaging methods quantitatively. A noise reduction of up to a factor of 5 for retardation and 3 for optic axis orientation could be achieved. Furthermore, the use of a temporal window for averaging Stokes vector elements instead of the previously used spatial window improved the spatial resolution of DOPU images by eliminating the convolution of the image structure with the evaluation window.
Interesting questions in the context of noise reduction by averaging are the contributions of statistical noise (shot noise, excess noise, receiver noise) and deterministic noise (speckle noise), the limit of the number of frames to be averaged, beyond which no further improvement is achieved, and the trade-off between S/N improvement by averaging and loss of spatial resolution. If the spatial relation between instrument and eye would be exactly the same for all averaged B-scans, each B-scan would have exactly the same speckle structure, and averaging would only reduce the statistic noise but not the deterministic speckle noise. To improve the speckle noise, averaging has to be performed over several independent speckles, i.e., over a volume that contains several speckles.
Different mechanisms can generate independent (decorrelated) speckles at the same nominal position (same settings of x-y scanner): (i) Rotations and axial motions of the eye. These rotations are translated into lateral shifts of the image, so these movements are equivalent to image shifts in three dimensions. (ii) A translational movement of the eye with respect to the illuminating beam will cause the beam to enter the pupil at different positions, which slightly changes the angle of the beam incident on the retina; oblique beams sample slightly different positions, giving rise to a different speckle structure. (iii) Focus changes by residual accommodation during the scan will change the wave fronts and might also slightly change the speckle pattern. The latter effect (iii) is probably negligible for an unaccommodated eye looking towards infinity, and within the short acquisition time of less than a second. Regarding point (ii): If we assume a maximum shift of the eye's entrance pupil of ~ ± 1 mm, and a focus position at the photoreceptor layer, the maximum lateral shift of the beam at the surface of a retina of ~300 µm thickness is of the order of ± 10 µm, i.e., in total not more than the size of a speckle (in transverse direction). This effect might therefore slightly add to averaging of independent speckles, however, we think that mechanism (i) is the dominating one.
Speckle reduction by averaging of independent speckles generated by mechanism (i) can only work if the motion correction is not perfect, i.e. if there is a shift between individual frames so that several independent speckles are covered by averaging. We can roughly estimate the number of independent speckles included by the averaging process. The shift of individual frames that the software motion correction algorithm applies can be extracted from the processing algorithm. It covers ~45-60 µm or ~3-4 speckles in x direction and ~10 µm (2 speckles) in z-direction. From the sharpness of individual very fine features (e.g. external limiting membrane, capillaries in anterior retinal layers), we can conclude that the motion correction algorithm reduces the residual motion artifacts to within ~0.5 – 1 speckle diameters (these residual artifacts are probably caused by in-frame distortions). Therefore, within the x-z plane ~2-4 independent speckles contribute to the averaging process. Previous work has shown, that the eye motions in y-direction have approximately similar magnitude as in x-direction [7], i.e. extends over ~3-4 speckles. Motions in y-direction are not corrected by our algorithm. Therefore, we estimate that a total of ~6 – 16 independent speckles are included in the averaging process. Since noise reduction scales with the square root of independent measurements, one would expect an improvement by roughly a factor of 3-4 from speckle noise reduction.
If we compare these results with Tables 1 and 2, we see that the total noise reduction achievable is somewhat larger, at least in the case of retardation imaging. For this quantity, the standard deviation can be improved by a factor of ~5 (averaging method 2) for averaging over n = 60 frames. The additional reduction can be explained by reduction of statistical noise. For the axis orientation, the improvement factor is only ~3. The reason for this discrepancy is probably that the axis orientation is encoded in the phase difference of the signals. Phase images suffer more from statistical noise than intensity images. A look at Table 3 shows that, while the maximum number of n beyond which no further noise reduction is observed, is n ~40 for retardation images, whereas axis imaging benefits from larger n, and with n = 60 the maximum number of useful frames might not yet be reached.
There is a trade-off between noise reduction by averaging and loss of spatial resolution. The averaging over several spatially separated speckles means that each resolved image point contains signal smeared out over the physical locations corresponding to these speckles, equivalent to a resolution loss. However, one should be aware that the definition of resolution within an OCT image of a scattering sample is not simply the size of the point spread function (as in case of imaging a single point scatterer or a specular reflecting surface). While the speckle size resembles the point spread function, a single speckle cannot be regarded as corresponding to a resolved structure. Usually a few speckles are necessary to judge, e.g., the position of a boundary. Therefore, the sampling over 2-4 independent speckles within the x-z plane seems not to be a real sacrifice, given the reduced noise. On the other hand, the averaging in y-direction over ~3-4 speckles is equivalent to acquiring data from a sheet of ~50-100 µm thickness, i.e., fine details perpendicular to the B-scan direction will get smeared out by the averaging procedure.
Finally, it should be mentioned that the numbers presented in this analysis are system dependent. Different system parameters (axial and transverse resolution, sampling density, imaging speed) can influence the motion correction and also the speckle size, and thereby the number of independent speckles that can be averaged. Furthermore, the software based motion correction works only with subjects of sufficient fixation capabilities. In the data sets used in this work, we observed translational shifts among all B-scans of a set within a range of 30 to 500µm, which were well tolerated by the motion correction algorithm. For larger motions that might occur with unstable patients, hardware tracking might be required.
Acknowledgment
Financial support from the Austrian Science Fund (FWF grant P19624-B02) and from the European Union project FUN OCT (FP7 HEALTH, contract no. 201880) is gratefully acknowledged.
References and links
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]
2. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]
3. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography (Springer, 2008).
4. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt. 7(3), 457–463 (2002). [CrossRef] [PubMed]
5. R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, “Performance of fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). [CrossRef] [PubMed]
6. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. Schuman, A. Kowalczyk, and J. S. Duker, “Three-dimensional retinal imaging with high-speed ultrahigh-resolution optical coherence tomography,” Ophthalmology 112(10), 1734–1746 (2005). [CrossRef] [PubMed]
7. B. Potsaid, I. Gorczynska, V. J. Srinivasan, Y. Chen, J. Jiang, A. Cable, and J. G. Fujimoto, “Ultrahigh speed spectral / Fourier domain OCT ophthalmic imaging at 70,000 to 312,500 axial scans per second,” Opt. Express 16(19), 15149–15169 (2008). [CrossRef] [PubMed]
8. M. R. Hee, D. Huang, E. A. Swanson, and J. G. Fujimoto, “Polarization-sensitive low-coherence reflectometer for birefringence characterization and ranging,” J. Opt. Soc. Am. B 9(6), 903–908 (1992). [CrossRef]
9. J. F. de Boer, T. E. Milner, M. J. C. van Gemert, and J. S. Nelson, “Two-dimensional birefringence imaging in biological tissue by polarization-sensitive optical coherence tomography,” Opt. Lett. 22(12), 934–936 (1997). [CrossRef] [PubMed]
10. J. F. de Boer, T. E. Milner, and J. S. Nelson, “Determination of the depth-resolved Stokes parameters of light backscattered from turbid media by use of polarization-sensitive optical coherence tomography,” Opt. Lett. 24(5), 300–302 (1999). [CrossRef] [PubMed]
11. S. Michels, M. Pircher, W. Geitzenauer, C. Simader, E. Götzinger, O. Findl, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Value of polarisation-sensitive optical coherence tomography in diseases affecting the retinal pigment epithelium,” Br. J. Ophthalmol. 92(2), 204–209 (2008). [CrossRef] [PubMed]
12. C. Ahlers, E. Götzinger, M. Pircher, I. Golbaz, F. Prager, C. Schütze, B. Baumann, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Imaging of the retinal pigment epithelium in age-related macular degeneration using polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 51(4), 2149–2157 (2010). [CrossRef] [PubMed]
13. B. Cense, T. C. Chen, B. H. Park, M. C. Pierce, and J. F. de Boer, “Thickness and birefringence of healthy retinal nerve fiber layer tissue measured with polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 45(8), 2606–2612 (2004). [CrossRef] [PubMed]
14. M. Yamanari, M. Miura, S. Makita, T. Yatagai, and Y. Yasuno, “Phase retardation measurement of retinal nerve fiber layer by polarization-sensitive spectral-domain optical coherence tomography and scanning laser polarimetry,” J. Biomed. Opt. 13(1), 014013 (2008). [CrossRef] [PubMed]
15. E. Götzinger, M. Pircher, B. Baumann, C. Hirn, C. Vass, and C. K. Hitzenberger, “Retinal nerve fiber layer birefringence evaluated with polarization sensitive spectral domain OCT and scanning laser polarimetry: a comparison,” J. Biophoton. 1(2), 129–139 (2008). [CrossRef] [PubMed]
16. M. Pircher, E. Götzinger, R. Leitgeb, H. Sattmann, O. Findl, and C. K. Hitzenberger, “Imaging of polarization properties of human retina in vivo with phase resolved transversal PS-OCT,” Opt. Express 12(24), 5940–5951 (2004). [CrossRef] [PubMed]
17. E. Götzinger, M. Pircher, and C. K. Hitzenberger, “High speed spectral domain polarization sensitive optical coherence tomography of the human retina,” Opt. Express 13(25), 10217–10229 (2005). [CrossRef] [PubMed]
18. M. Pircher, E. Götzinger, O. Findl, S. Michels, W. Geitzenauer, C. Leydolt, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Human macula investigated in vivo with polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 47(12), 5487–5494 (2006). [CrossRef] [PubMed]
19. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography,” Opt. Express 16(21), 16410–16422 (2008). [CrossRef] [PubMed]
20. M. Miura, M. Yamanari, T. Iwasaki, A. E. Elsner, S. Makita, T. Yatagai, and Y. Yasuno, “Imaging polarimetry in age-related macular degeneration,” Invest. Ophthalmol. Vis. Sci. 49(6), 2661–2667 (2008). [CrossRef] [PubMed]
21. B. Cense, W. H. Gao, J. M. Brown, S. M. Jones, R. S. Jonnal, M. Mujat, B. H. Park, J. F. de Boer, and D. T. Miller, “Retinal imaging with polarization-sensitive optical coherence tomography and adaptive optics,” Opt. Express 17(24), 21634–21651 (2009). [CrossRef] [PubMed]
22. C. E. Saxer, J. F. de Boer, B. H. Park, Y. Zhao, Z. Chen, and J. S. Nelson, “High-speed fiber based polarization-sensitive optical coherence tomography of in vivo human skin,” Opt. Lett. 25(18), 1355–1357 (2000). [CrossRef] [PubMed]
23. M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express 16(8), 5892–5906 (2008). [CrossRef] [PubMed]
24. W. Y. Oh, S. H. Yun, B. J. Vakoc, M. Shishkov, A. E. Desjardins, B. H. Park, J. F. de Boer, G. J. Tearney, and B. E. Bouma, “High-speed polarization sensitive optical frequency domain imaging with frequency multiplexing,” Opt. Express 16(2), 1096–1103 (2008). [CrossRef] [PubMed]
25. M. K. Al-Qaisi and T. Akkin, “Polarization-sensitive optical coherence tomography based on polarization-maintaining fibers and frequency multiplexing,” Opt. Express 16(17), 13032–13041 (2008). [CrossRef] [PubMed]
26. E. Götzinger, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Polarization maintaining fiber based ultra-high resolution spectral domain polarization sensitive optical coherence tomography,” Opt. Express 17(25), 22704–22717 (2009). [CrossRef] [PubMed]
27. H. Wang, M. K. Al-Qaisi, and T. Akkin, “Polarization-maintaining fiber based polarization-sensitive optical coherence tomography in spectral domain,” Opt. Lett. 35(2), 154–156 (2010). [CrossRef] [PubMed]
28. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]
29. N. Iftimia, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography by “path length encoded” angular compounding,” J. Biomed. Opt. 8(2), 260–263 (2003). [CrossRef] [PubMed]
30. M. Pircher, E. Gotzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8(3), 565–569 (2003). [CrossRef] [PubMed]
31. A. E. Desjardins, B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Speckle reduction in OCT using massively-parallel detection and frequency-domain ranging,” Opt. Express 14(11), 4736–4745 (2006). [CrossRef] [PubMed]
32. A. Ozcan, A. Bilenca, A. E. Desjardins, B. E. Bouma, and G. J. Tearney, “Speckle reduction in optical coherence tomography images using digital filtering,” J. Opt. Soc. Am. A 24(7), 1901–1910 (2007). [CrossRef] [PubMed]
33. A. Wong, A. Mishra, K. Bizheva, and D. A. Clausi, “General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery,” Opt. Express 18(8), 8338–8352 (2010). [CrossRef] [PubMed]
34. U. E. K. Wolf-Schnurrbusch, L. Ceklic, C. K. Brinkmann, M. E. Iliev, M. Frey, S. P. Rothenbuehler, V. Enzmann, and S. Wolf, “Macular thickness measurements in healthy eyes using six different optical coherence tomography instruments,” Invest. Ophthalmol. Vis. Sci. 50(7), 3432–3437 (2009). [CrossRef] [PubMed]
35. C. K. Hitzenberger, E. Goetzinger, M. Sticker, M. Pircher, and A. F. Fercher, “Measurement and imaging of birefringence and optic axis orientation by phase resolved polarization sensitive optical coherence tomography,” Opt. Express 9(13), 780–790 (2001). [CrossRef] [PubMed]
36. M. Pircher, E. Götzinger, H. Sattmann, R. A. Leitgeb, and C. K. Hitzenberger, “In vivo investigation of human cone photoreceptors with SLO/OCT in combination with 3D motion correction on a cellular level,” Opt. Express 18(13), 13935–13944 (2010). [CrossRef] [PubMed]
37. A. Szkulmowska, M. Szkulmowski, A. Kowalczyk, and M. Wojtkowski, “Phase-resolved Doppler optical coherence tomography-limitaions and improvemens,” Opt. Express 33, 1425–1427 (2008).
38. V. X. D. Yang, M. L. Gordon, B. Qi, J. Pekar, S. Lo, E. Seng-Yue, A. Mok, B. C. Wilson, and I. Vitkin, “High speed, wide velocity dynamic range Doppler optical coherence tomography (part I): system design, signal processing, and performance,” Opt. Express 11(7), 794–809 (2003). [CrossRef] [PubMed]
39. B. Fortune, G. A. Cull, and C. F. Burgoyne, “Relative course of nerve fiber layer birefringence and thickness and retinal function changes after optic nerve head transsection,” Invest. Ophthalmol. Vis. Sci. 49(10), 4444–4452 (2008). [CrossRef] [PubMed]
40. M. Pircher, E. Götzinger, B. Baumann, and C. K. Hitzenberger, “Corneal birefringence compensation for polarization sensitive optical coherence tomography of the human retina,” J. Biomed. Opt. 12(4), 041210 (2007). [CrossRef] [PubMed]
41. K. Wiesauer, M. Pircher, E. Goetzinger, C. K. Hitzenberger, R. Engelke, G. Ahrens, G. Gruetzner, and D. Stifter, “Transversal ultrahigh-resolution polarizationsensitive optical coherence tomography for strain mapping in materials,” Opt. Express 14(13), 5945–5953 (2006). [CrossRef] [PubMed]
42. S. Jiao, R. Knighton, X. Huang, G. Gregori, and C. Puliafito, “Simultaneous acquisition of sectional and fundus ophthalmic images with spectral-domain optical coherence tomography,” Opt. Express 13(2), 444–452 (2005). [CrossRef] [PubMed]
43. E. Götzinger, M. Pircher, B. Baumann, C. Hirn, C. Vass, and C. K. Hitzenberger, “Analysis of the origin of atypical scanning laser polarimetry patterns by polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 49(12), 5366–5372 (2008). [CrossRef] [PubMed]
44. A. W. Dreher, K. Reiter, and R. N. Weinreb, “Spatially resolved birefringence of the retinal nerve fiber layer assessed with a retinal laser ellipsometer,” Appl. Opt. 31(19), 3730–3735 (1992). [CrossRef] [PubMed]
