## Abstract

We propose a transit-time based method to ascertain the azimuth angle of a velocity vector by spectral-domain Doppler optical coherence tomography (DOCT), so that three-dimensional (3-D) velocity vector can be quantified. A custom-designed slit plate with predetermined slit orientation is placed into the sample beam to create three delay-encoded sub-beams of different beam shape for sample probing. Based on the transit-time analysis for Doppler bandwidth, the azimuth angle within 90° range is evaluated by exploitation of the complex signals corresponding to three path length delays. 3-D velocity vector is quantified through further estimating of Doppler angle and flow velocity by combined Doppler shift and Doppler bandwidth measurements. The feasibility of the method is demonstrated by good agreement between the determined azimuth angles and the preset ones, and further confirmed by velocity vector measurement of flowing solution inside a capillary tube.

©2010 Optical Society of America

## 1. Introduction

Doppler optical coherence tomography (DOCT) is a powerful technique for determining the localized flow velocity in highly scattering media such as biological tissues [1]. Axial component of the velocity vector is directly detectable in a DOCT, but the other two components of the 3-D velocity vector have been ascertained through ingenious approaches. Dave and Milner [2] were the first to measure Doppler angle by detecting velocity at two different angles with orthogonal polarization modes, where two Wollaston prisms were adopted with a prerequisite of polarization-insensitive property of the sample. Pedersen et al [3] used a glass plate to create delay-encoded sub-beams with two incident angles for determining Doppler angle. Doppler spectrum broadening [4] was exploited to measure transversal velocity and hence the Doppler angle. However, the azimuth angle of the velocity vector was still undetermined. Full determination of velocity vector was realized by a four-channel quadrant detector in a free-space system [5], but such implementation is not applicable to a fiber-based system. A beam divider was proposed to quantify the 3-D velocity vector in a fiber-based spectral-domain DOCT system [6]. The beam divider encoded the probe beam with five different path length delays requiring too much depth range, which causes severe restriction on spectral-domain OCT system. Other groups [7–9] tried to get the velocity vector from the estimation of the flowing direction based on reconstructed 3-D orientation of a blood vessel. However, the velocity direction of different position inside the vessel is not strictly consistent with the vessel orientation.

In this paper, we present a transit-time based method in a fiber-based spectral-domain DOCT system to ascertain the azimuth angle of the velocity vector. The method is especially useful to determine the azimuth angle with high accuracy within 90° range under prior knowledge of azimuth angle known from the reconstructed 3-D orientation of the blood vessel. By inserting a slit plate with predetermined slit orientation into the sample arm, the probe beam is divided into three delay-encoded sub-beams with different beam shape. Three independent complex signals corresponding to different path length delays are obtained by Fourier transform of the detected spectral interferogram. The azimuth angle within 90° range is then determined by Doppler bandwidth measured from the three complex signals. 3-D velocity vector is thus quantified by further estimating of Doppler angle and flow velocity by combined Doppler shift and Doppler bandwidth measurements [10]. The available depth is 1/3 of the maximum imaging depth in the proposed method, in contrast to 1/5 of the maximum imaging depth in the previously reported method [6]. Therefore, the demand on imaging depth for spectral domain DOCT system is relaxed. Higher accuracy on determining velocity vector is hence envisioned due to the depth dependent sensitivity fall-off in spectral-domain OCT. Finally, feasibility of the proposed method is demonstrated by phantom experiments.

## 2. Theory

The method is based on the full exploitation of Doppler spectrum broadening. There are a number of sources, including flow velocity, turbulence, Brownian motion, as well as the probe beam geometry, contribute to the broadened Doppler spectrum [4]. Doppler spectrum broadening is mainly determined by the transit time of scattering particles passing through the focused probe beam [11]. The transit time in the focal area is

Doppler bandwidth*B*in evaluation of the Doppler spectrum broadening is the inverse of the transit time, i.e.,As illustrated in Fig. 1 , V is the flow velocity,

*α*is the Doppler angle between the flowing direction and the optical axis of the probe beam,

*w*denotes the effective waist diameter of the probe beam in the focal area, which is the transit length of moving particles.

*B*

_{0}accounts for the contributions from other sources that are independent of the macroscopic flow velocity.

To quantify the 3-D velocity vector of moving particles, it is required to obtain three quantities, which are flow velocity *V*, Doppler angle *α* and azimuth angle *φ*. As shown in Fig. 1(a), the slit plate is made of glass with air slit in the middle. The glass part of the plate has a thickness of *t* with refractive index of *n _{g}*. The slit with a width of

*d*is oriented along x axis. The plate is introduced in the sample arm of a spectral-domain DOCT system between the collimator and galvo mirror shown in Fig. 1(b). Collimated beam with a diameter of

*D*illuminates the plate and focused to the scattering particles by a focusing lens with the focal length of

*f*. Backscattered light from the moving particles is collected by the focusing lens and passes through the plate again. In double passing of the plate, the light could go through the plate via air slit first and then return via air slit or glass. Alternatively, the light could go through the plate via glass first and then return via air slit or glass. Therefore, the slit plate makes four different beam paths: air to air (AA), air to glass (AG), glass to air (GA) and glass to glass (GG). Since AG and GA have the same optical path delay, there are three independent path length delays, which can be denoted by 0,

*t*(

*n*-1) and 2

_{g}*t*(

*n*-1). In the following text, the double-subscript of the transit length

_{g}*w*and the Doppler bandwidth

*B*indicates the corresponding beam path.

The maximum imaging depth in optical length of the spectrometer based spectral-domain DOCT system is determined by *Mλ*
^{2}/4Δ*λ*, where *M* is the number of pixels of the CCD, *λ* is the central wavelength of the light source, and Δ*λ* is the full spectral range of the spectrometer. In the previously reported method [6], neighbouring sub-images of the sample with an optical depth larger than *Mλ*
^{2}/20Δ*λ* will overlap, while in the proposed method this threshold is improved to be *Mλ*
^{2}/12Δ*λ*. Therefore, the demand on imaging depth for the spectral-domain DOCT system is relaxed.

Figure 2(a) and (b)
shows the simulation results of normalized intensity patterns in the focal plane of the focusing lens corresponding to collimated illumination via air slit (*d*=1 mm) and glass, respectively. *x* and *y* are horizontal and vertical coordinates of the focal plane. To analyze the transit time, the beam size is estimated from the contour of the central maximum of the intensity pattern. Under illumination via air slit demonstrated in Fig. 2(a), the width of the central maximum along *x* direction and *y* direction are 2*λf* /*D* and 2*λf* /*d*, respectively. The fitted contour of the central maximum is almost a rectangle. Therefore, the transit length *w*
_{AA} that a particle passes through the central pattern with an azimuth angle *φ* can be expressed as

*φ*approaches 90°,

*w*

_{AA}becomes very large, a small Doppler bandwidth comparable to the background

*B*

_{0}is resulted and the information on azimuth angle

*φ*is impossible to recover accurately.

Under illumination via glass as demonstrated in Fig. 2(b), contour of the central maximum can be fitted by elliptic equation. The fitted contours with different width of *d* are plotted in Fig. 2(c). Obviously, the central maximum under illumination via glass approaches an Airy pattern as the air slit *d* decreases. The transit length *w*
_{GG} that a particle passes through this central pattern with an azimuth angle *φ* is given by 2(*x*
_{0}
^{2}cos^{2}
*φ*+*y*
_{0}
^{2}sin^{2}
*φ*)^{1/2}, where *x*
_{0} and *y*
_{0} represent the semimajor axis and semiminor axis of the fitted ellipse, respectively. For instance, with system parameters of *λ*=0.84 μm, *f*=75 mm, *D*=5.2 mm and *d*=1 mm, the fitted elliptic contour of the central maximum is (x/15.77)^{2}+(y/11.50)^{2}=1. Thus we have *x*
_{0}=15.77 μm and *y*
_{0}=11.50 μm.

The azimuth angle within 90° range is ascertained by full exploitation of complex signals corresponding to three path length delays. In spectral-domain DOCT, a depth profile is obtained by calculating the Fourier transform of the spectral interferogram in the k-space. The depth profile can be recognized as a complex signal represented by $\tilde{\Gamma}(z)=A(z)\mathrm{exp}[i\Phi (z)]$, where *z* is the axial coordinate, *A*(*z*) and *Φ*(*z*) denote the amplitude and phase terms, respectively. *A*(*z*) is used for structural information reconstruction, and phase *Φ*(*z*) is used for Doppler information construction. Then $\tilde{\Gamma}(z)$ is separated into three complex signals corresponding to three path length delays. Assuming ${\tilde{\Gamma}}_{AA}(z)$ is the complex signal corresponding to sub-beam AA, ${\tilde{\Gamma}}_{AA}(z)$ is extracted from $\tilde{\Gamma}(z)$ with z in the range of (z_{0}-*δ*/2, z_{0}+*δ*/2), where *δ*=*t*(*n _{g}*-1)/2 and z

_{0}represents the axial position corresponding to the center of sample depth. ${\tilde{\Gamma}}_{AG+GA}(z+\delta )$ is the complex signal corresponding to sub-beams AG and GA, which adds a delay shift of

*δ*in axial direction compared with ${\tilde{\Gamma}}_{AA}(z)$. Similarly, ${\tilde{\Gamma}}_{GG}(z+2\delta )$ is the signal corresponding to sub-beam GG, which adds a delay shift of 2

*δ*in axial direction compared with ${\tilde{\Gamma}}_{AA}(z)$.

By coherently summation of the complex signals with three path length delays, we have

*B*

_{sum}, is independent on the azimuth angle

*φ*because of circular symmetry of the equivalent clear aperture. The transit length corresponding to the summed complex signal,

*w*

_{sum}, is the diameter of well-known Airy pattern given by 8

*λf*/π

*D*irrespective of the azimuth angle.

According to Eq. (2), when the background *B*
_{0} is subtracted, the ratio of the Doppler bandwidth obtained from ${\tilde{\Gamma}}_{AA}(z)$ to that obtained from ${\tilde{\Gamma}}_{sum}(z)$ is given by

Similarly, the ratio of the Doppler bandwidth obtained from ${\tilde{\Gamma}}_{GG}(z+2\delta )$ to that obtained from ${\tilde{\Gamma}}_{sum}(z)$ is given by

Based on the relationship between Doppler bandwidth *B* and standard deviation *σ*, i.e. $B=4\sigma $, the Doppler bandwidth is obtained from the standard deviation, which could be calculated from complex signals of successive A-scans according to [4]:

*T*is the time interval between two A-scans,

*N*is the number of A-lines used for averaging, ${\tilde{\Gamma}}_{j}(z)$ and ${\stackrel{~}{\Gamma}}_{j+1}(z)$ are the complex signals corresponding to the

*j*th A-scan and its next A-scan, ${\stackrel{~}{\Gamma}}_{j}^{*}(z)$ is the conjugate of ${\tilde{\Gamma}}_{j}(z)$.

After calculating Doppler bandwidth *B*
_{AA}, *B*
_{GG} and *B*
_{sum}, azimuth angle *φ* within 90° range can be determined using Eq. (5) or (6). However, the determination of an unknown *φ* through Eq. (5) is under a prior knowledge of its range, which can be obtained by Eq. (6). The other two quantities required for the quantification of the velocity vector, Doppler angle *α* and flow velocity *V*, are estimated by combined Doppler shift and Doppler bandwidth measurements from the coherently summed complex signal ${\tilde{\Gamma}}_{sum}(z)$ based on previously reported method [10].

The relationship between Doppler shift *f _{d}* and flow velocity

*V*is

*λ*is the central wavelength in the media.

Doppler shift *f _{d}* is calculated by phase difference of successive A-scans:

From Eq. (2) and Eq. (8), the Doppler angle *α* is calculated by

*V*is obtained by

The flow chart of signal processing procedure is shown in Fig. 3 .

Intuitively, the case via air slit shown in Fig. 2(a) is more sensitive to slit orientation in comparison with that via glass shown in Fig. 2(b), and higher accuracy on determination of azimuth angle might be envisioned using Eq. (5) than using Eq. (6). However, the SNR corresponding to ${\tilde{\Gamma}}_{AA}(z)$ is worse than ${\tilde{\Gamma}}_{GG}(z+2\delta )$ due to the small width of air slit *d*. Therefore, both Eq. (5) and (6) are tried for the evaluation of azimuth angle in the following experiments.

## 3. Experiments

To validate the feasibility of the transit-time based method for 3-D velocity vector quantification, the fiber-based spectral-domain OCT system previously described [13, 14] is functionally extended to a DOCT system besides inserting the slit plate into the sample arm as shown in Fig. 1(b). The light source is a superluminescent diode with a central wavelength of 840 nm and bandwidth of 45 nm. In the sample arm, the 5.2 mm diameter beam passes through the slit plate, scanned by a galvo mirror and focused by a 75 mm focal length lens. In the spectrometer, the interfered beam is collimated by a 60 mm focal length lens, dispersed by a volume phase holographic grating of 1200 lines/mm, focused by a 150 mm focal length lens, and detected by a line-CCD camera of 2048 pixels with 20 KHz A-scan rate. Using a 1.2 mW sample beam, the sensitivity measured near DC is 108 dB with an axial resolution of 7.5 µm in air. The dynamic range is 68 dB. There is a 16 dB falloff in sensitivity over 2 mm depth. The unambiguous velocity range lies on the maximum detectable Doppler shift, which is restricted by the A-scan rate. The determinable Doppler angle is from 0° to 180°. The determinable azimuth angle is within 90° range inside a quadrant.

0.25% aqueous solution of polystyrene beads (0.53 μm diameter) in a capillary tube with inner diameter of 150 μm is used as the sample. Flowing speed of the scattering beads within the tube is controlled by a syringe pump. The tube orientation is set at azimuth angle *φ* starting from *x* axis and Doppler angle *α* starting from *z* axis as depicted in Fig. 1. The air silt width *d* is set to be 1 mm or 0.8 mm for comparable signal intensities of three sub-images corresponding to the three path length delays. The glass part of the plate is made of BK7, which has a thickness of 1.7 mm with a refractive index of 1.51 at 830 nm wavelength. In our experiment, group velocity mismatch among the different paths introduced by the slit plate is not considered since the thickness of the glass part is relatively small.

Figure 4 demonstrates the cross-sectional Doppler images representing Doppler shift and Doppler bandwidth of the capillary tube with flowing solution of polystyrene beads inside. The three sub-images of the same sample corresponding to three path length delays are obtained simultaneously by the fiber-based spectral-domain DOCT system.

To investigate the ability of the proposed method to determine azimuth angle, experiments corresponding to ten different preset azimuth angles within 90° range are conducted with slit width *d* of 1 mm and 0.8 mm. Under the condition that the orientation of capillary tube at a fixed Doppler angle of 76°, the volume velocity of flowing solution is varied from 1 μl/min to 5 μl/min at increment step of 1 μl/min, equivalent to an average velocity changing from 0.94 mm/s to 4.72 mm/s.

M-mode data acquisition is to perform several A-scans at a transversal position without transversal scanning [4]. In the experiment, the probe beam is focused to the center of the capillary tube. *N* A-scans are performed for averaging (*N*=1000). With A-scan rate at 20 KHz, the time needed for one M-mode data acquisition is 50 ms. The reason why we use M-mode in the experiment is to exclude other potential affecting factors which include: 1) deflexion of incident ray caused by geometric shape of capillary tube and the refraction index mismatch; 2) beam shape variation with the position of the scanning beam due to optical aberrations; 3) additional phase noise induced by the transversal scanning.

After calculating the Doppler bandwidth obtained from ${\tilde{\Gamma}}_{AA}(z)$, ${\tilde{\Gamma}}_{GG}(z+2\delta )$ and ${\tilde{\Gamma}}_{sum}(z)$ under different flow velocities, the corresponding background *B*
_{0} can be deduced through linear fitting and extending to zero velocity. The estimated background *B*
_{0} for cases of ${\tilde{\Gamma}}_{AA}(z)$, ${\tilde{\Gamma}}_{GG}(z+2\delta )$ and ${\tilde{\Gamma}}_{sum}(z)$ with *d*=1 mm are 128.8, 140.4 and 136.0 Hz, respectively, and with *d*=0.8 mm, those are 121.6, 139.2 and 133.2 Hz, respectively.

The azimuth angles estimated by Eq. (5) and Eq. (6) with preset flow velocity of 4.72 mm/s are shown in Fig. 5(a) and (b)
, respectively. As shown in Fig. 5(a), the estimated azimuth angles are in good agreement with the preset ones for *φ* smaller than 80°. The accuracy for two values of *d* is comparable. Invalid estimation of azimuth angle occurs when it approaches 90° due to the big influence of Doppler bandwidth background, thus the calculated data at 80° and 90° are discarded. Compared with the results shown in Fig. 5(a), the estimated azimuth angles in Fig. 5(b) are in better agreement with the preset ones especially for *φ* near 90°. The data for *φ* around 45° is in perfect match with the preset ones, since *w*
_{GG} as a function of *φ* is approximately linear in this situation. However, the data for *φ* near 0° or 90° is not accurate for the reason that *w*
_{GG} is insensitive to the variation of *φ* in these two situations. As can be seen, the result for *d*=1 mm is in better agreement with the preset ones compared with *d*=0.8 mm. As *d* becomes larger, the fitted elliptic contour of the central maximum under illumination via glass becomes more elongated (see Fig. 2(c)). Therefore the transit length *w*
_{GG} is more sensitive to the variation of *φ*. The data of 0° is discarded since the calculated cos^{2}
*φ* (see Eq. (6)) is invalid owing to added noise. The RMS errors of the estimated angles corresponding to Fig. 5(a) and (b) for two values of *d* are summarized in Table 1
. The difference between estimated azimuth angles based on Eq. (5) and preset ones is within the range of 0.29° to 9.37°, while based on Eq. (6) the difference is within the range of 0.19° to 12.90°. Higher accuracy of estimation on azimuth angle is achieved using Eq. (5) than using Eq. (6) if invalid estimations at 80° and 90° are excluded. On the other hand, method based on Eq. (6) can still obtain effective estimation of azimuth angle even it approaches 90° due to relatively high SNR.

After azimuth angle *φ* is determined, the Doppler shift and Doppler bandwidth obtained from the coherently summed complex signal ${\tilde{\Gamma}}_{sum}(z)$ are used for the estimation of flow velocity *V* and Doppler angle *α*. The estimated results are shown in Fig. 6(a) and (b)
, where the capillary tube is positioned at azimuth angle of 45°. It can be seen that the estimated data for *d*=1 mm and 0.8 mm are both in good agreement with the preset ones under different flow velocities. The estimated flow velocity vector is thus quantified and plotted in Fig. 7
with the preset flow velocity of 4.72 mm/s.

In the text above, the Doppler bandwidths are computed for in focus situation. In the situation out of focus, the actual transit length *w* is enlarged due to bigger lateral spot size. Therefore, according to the relationship of Doppler bandwidth and transit length in Eq. (2), the Doppler bandwidth *B* decreases. In our method, the azimuth angle *φ* is obtained by calculating the ratio of *B*
_{AA} and *B*
_{sum} or *B*
_{GG} and *B*
_{sum} according to Eq. (5) and (6). Although the Doppler bandwidth decreases, the variation of *B*
_{AA}/*B*
_{sum} or *B*
_{GG}/*B*
_{sum} is usually negligible since the decreasing of *B*
_{AA}, *B*
_{GG} and *B*
_{sum} is almost proportional, especially when the defocusing distance is small.

However, defocusing effect will result in an incorrect estimation of Doppler angle. The Doppler angle is calculated by the combination of Doppler bandwidth and Doppler shift as demonstrated by Eq. (10). Where *w* is the theoretical transit length, *f _{d}* is the Doppler shift which is not affected by defocusing. Therefore, the Doppler angle is underestimated.

In the sample arm of our system, the 5.2 mm diameter beam is focused by a 75 mm focal length lens, thus the depth of focus is 445 μm. The depth of focus is much larger than the optical depth of the experimental phantom. The probe beam is weakly focused, and in this case, the defocusing effect is controlled to the minimum.

## 4. Conclusion

In summary, we present a method to quantify a velocity vector in spectral-domain DOCT based on transit-time analysis of scattering particles passing through the focusing pattern formed by delay-encoded sub-beams. Doppler bandwidth measured from the complex signals corresponding to three different delays as well as the coherently summed complex signal are fully exploited for the determination of azimuth angle. From conducted experiments on azimuth angle estimation of a capillary tube at preset orientations, we compare the accuracy of two methods based on Doppler bandwidth ratio measured from different complex signals. Velocity vector of flowing solution inside a capillary tube quantified with discernible velocity profile is in good agreement with expectation. Under prior knowledge of azimuth angle known from the reconstructed 3-D orientation of the vessel, it is feasible to determine the azimuth angle with high accuracy within 90° range and thus obtain velocity profile within the vessel. Further study on the numerical dispersion compensation is necessary in the application of highly dispersive biomedical sample.

## Acknowledgement

This work was supported by National High Technology Research and Development Program of China (2006AA02Z4E0) and Natural Science Foundation of China (60878057, 60978037).

## References and links

**1. **Z. Chen, Y. Zhao, S. M. Srinivas, J. S. Nelson, N. Prakash, and R. D. Frostig, “Optical Doppler Tomography,” IEEE J. Sel. Top. Quantum Electron. **5**(4), 1134–1142 (1999). [CrossRef]

**2. **D. P. Davé and T. E. Milner, “Doppler-angle measurement in highly scattering media,” Opt. Lett. **25**(20), 1523–1525 (2000). [CrossRef]

**3. **C. J. Pedersen, D. Huang, M. A. Shure, and A. M. Rollins, “Measurement of absolute flow velocity vector using dual-angle, delay-encoded Doppler optical coherence tomography,” Opt. Lett. **32**(5), 506–508 (2007). [CrossRef] [PubMed]

**4. **H. Ren, K. M. Brecke, Z. Ding, Y. Zhao, J. S. Nelson, and Z. Chen, “Imaging and quantifying transverse flow velocity with the Doppler bandwidth in a phase-resolved functional optical coherence tomography,” Opt. Lett. **27**(6), 409–411 (2002). [CrossRef]

**5. **A. Royset, T. Storen, F. Stabo-Eeg, and T. Lindmo, “Quantitative measurements of flow velocity and direction using Transversal Doppler Optical Coherence Tomography,” Proc. SPIE **6079**, 607925 (2006). [CrossRef]

**6. **Y. C. Ahn, W. Jung, and Z. Chen, “Quantification of a three-dimensional velocity vector using spectral-domain Doppler optical coherence tomography,” Opt. Lett. **32**(11), 1587–1589 (2007). [CrossRef] [PubMed]

**7. **R. Michaely, A. H. Bachmann, M. L. Villiger, C. Blatter, T. Lasser, and R. A. Leitgeb, “Vectorial reconstruction of retinal blood flow in three dimensions measured with high resolution resonant Doppler Fourier domain optical coherence tomography,” J. Biomed. Opt. **12**(4), 041213 (2007). [CrossRef] [PubMed]

**8. **Y. Wang, B. A. Bower, J. A. Izatt, O. Tan, and D. Huang, “In vivo total retinal blood flow measurement by Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. **12**(4), 041215 (2007). [CrossRef] [PubMed]

**9. **S. Makita, T. Fabritius, and Y. Yasuno, “Quantitative retinal-blood flow measurement with three-dimensional vessel geometry determination using ultrahigh-resolution Doppler optical coherence angiography,” Opt. Lett. **33**(8), 836–838 (2008). [CrossRef] [PubMed]

**10. **D. Piao, L. L. Otis, and Q. Zhu, “Doppler angle and flow velocity mapping by combined Doppler shift and Doppler bandwidth measurements in optical Doppler tomography,” Opt. Lett. **28**(13), 1120–1122 (2003). [CrossRef] [PubMed]

**11. **V. L. Newhouse, E. S. Furgason, G. F. Johnson, and D. A. Wolf, “The dependence of ultrasound Doppler bandwidth on beam geometry,” IEEE Trans. Sonics Ultrason. **SU-27**, 50–59 (1980).

**12. **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]

**13. **K. Wang, Z. Ding, T. Wu, C. Wang, J. Meng, M. Chen, and L. Xu, “Development of a non-uniform discrete Fourier transform based high speed spectral domain optical coherence tomography system,” Opt. Express **17**(14), 12121–12131 (2009), http://www.opticsinfobase.org/abstract.cfm?uri=oe-17-14-12121. [CrossRef] [PubMed]

**14. **K. Wang, Z. Ding, Y. Zeng, J. Meng, and M. Chen, “Sinusoidal B-M method based spectral domain optical coherence tomography for the elimination of complex-conjugate artifact,” Opt. Express **17**(19), 16820–16833 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-17-19-16820. [CrossRef] [PubMed]