## Abstract

For Fourier domain optical coherence tomography any sample movement during camera integration causes blurring of interference fringes and as such reduction of sensitivity for flow detection. The proposed method overcomes this problem by phase-matching a reference signal to the sample motion. The interference fringes corresponding to flow signal will appear frozen across the detector whereas those of static sample structures will be blurred resulting in enhanced contrast for blood vessels. An electro-optic phase modulator in the reference arm, driven with specific phase cycles locked to the detection frequency, allows not only for qualitative but also for quantitative flow detection already from the relative signal intensities. First applications to extract in-vivo retinal flow and to visualize 3D vascularization, i.e. optical vivisection, are presented.

©2007 Optical Society of America

## 1. Introduction

Fourier or spectral domain optical coherence tomography (FDOCT) has nowadays largely replaced time domain OCT due to its clear sensitivity and acquisition speed advantage [1-3]. This is due to the fact that the full sample depth structure is obtained in a parallel way: it is encoded as depth dependent modulation frequency of the light source spectrum. For the depth profile (A-scan) reconstruction it is in principle sufficient to perform a Fourier analysis of the spectral interference pattern that can be recorded in a single camera shot [4-6]. The speed advantage is apparent since no mechanical depth scanning is needed. In particular for fast retinal in-vivo imaging with high resolution in 3D, FDOCT has become the method of choice [7-9]. Recent developments enhance the clinical and biomedical potential of FDOCT by aiming from purely structural to functional tissue imaging revealing tissue dynamics and physiology. The imaging space is becoming highly multi-dimensional including polarization contrast [10-13], Doppler flow [13-17], and spectroscopic contrast [18,19].

Several studies of retinal blood flow using Laser Doppler flowmetry (LDF) outlined the importance of measuring vessel flow properties since it is an early indicator of retinal pathologies like glaucoma, diabetic retinopathy or age related macula degeneration [20-24]. Color Doppler FDOCT based on phase sensitive signal detection has been already applied to visualize retinal perfusion with high velocity sensitivity and to contrast for example systolic and diastolic blood flow [15]. The advantage of Doppler OCT as compared to LDF is the ability to localize retinal perfusion in depth [25]. Recent work applied elaborate image processing techniques to extract the retinal vascularization and to differentiate between choroidal and retinal flow bed [17]. The standard method to visualize retinal and choroidal perfusion is still fluorescence angiography. Such methods are invasive since in general fluorescein or indocyanine green (ICG) need to be injected to produce the contrasting of vascular structure on the ocular fundus. Although the adverse effects of such contrast agents could be reduced, a completely non-invasive assessment would still increase the comfort of the patient together with the easiness in handling such instruments. In addition, a fast tomographic method such as FDOCT offers the possibility to completely assess the retinal vascularization in 3D.

Still, spectrometer-based FDOCT suffers from inherent imaging artifacts such as complex ambiguity, aliasing, signal decay in depth, and last but not least interference fringe blurring. A recent approach to reduce the latter effect for retinal imaging has been the application of pulsed illumination [26]. Fringe blurring in particular limits the detection of moving structures such as flow. Figure 1 demonstrates a typical high resolution retinal tomogram with visible blood vessels. One can clearly observe vessels with signal from internal blood flow as well as vessels that appear empty. Fast flows occur especially in the region of the optic nerve head (ONH) where the central retinal artery and vein enter and leave the eye bulb. Hence, it is of great interest to keep the flow signal in this particular region for assessing retinal perfusion.

The present approach circumvents the limitation of interference fringe blurring by phase-matching the detection to the sample movement, similar to a photographer that moves the camera with the object. In effect, static structures will be blurred whereas moving parts such as flow will be enhanced. The present article gives first a detailed discussion of the theory of fringe blurring in FDOCT using broad bandwidth light sources. Subsequently the theory of resonant Doppler imaging will be introduced together with the realization of employing electro-optic phase modulation. It will be shown how flow can be evaluated not only qualitatively but also quantitatively on a pure intensity basis. It should be mentioned that an alternative method of phase independent Doppler OCT has been introduced recently by analysis of speckle patterns [27]. Finally the method is applied to measure retinal blood flow and 3D vascularization in the ONH region.

## 2. Theory

#### 2.1. Chromatic interference fringe blurring

It is a well known fact that in spectrometer-based FDOCT any sample movement causes interference fringe blurring and in consequence reduction of the signal-to-noise ratio (SNR). First analytic discussion of this effect has been given [28], but without accounting for the intrinsic chromaticity of FDOCT and without giving an analytic expression for the true FDOCT signal after Fourier transform. As we will see in due course, the fringe blurring causes not only loss of signal contrast but also signal peak broadening.

Let us assume, without loss of generalization, that we have a single reflecting interface in the sample arm. The reference field *E _{R}(k)* and the sample field

*E*are expressed as functions of the wavenumber

_{S}(k)*k*. We further assume, in the most general case, a time dependent phase shift

*ϕ(t)*between the sample and the reference field. The recorded spectral interference signal, written as number of generated photo-electrons, becomes

with *I _{S,R}*=

*E*

_{S,R}(k)E_{S,R}(k)^{*}, the asterisk denoting complex conjugate,

*z*

_{0}being the optical path length difference between sample and reference interface at

*t*=0,

*τ*the camera exposure time, and the proportionality factor

*γ*accounting for the detector conversion efficiency. Subsequently we will only consider the last part of Eq. (1) that actually contains information about the relative sample position. Denoting this part

*N*and assuming a phase shift

_{AC}(k)*ϕ(t)*that is linear in time, the integration yields [28]

$${N}_{\mathit{AC}}\left(k,\mathrm{\Delta \Phi}=0\right)=2\mathit{\gamma \tau}\sqrt{{I}_{S}\left(k\right)}\sqrt{{I}_{R}\left(k\right)}\mathrm{cos}\left(2k{z}_{0}\right),$$

where ΔΦ=*ϕ(τ)* is the total phase change during the integration time, and sinc(*x*)≡sin(*πx*)/(*πx*). The varying phase *ϕ* reduces the modulation depth of the spectral interference signal across the detector array. The origin of the phase shift can be manifold: often we are faced with involuntary sample movements, in particular during in-vivo measurements, or mechanical noise that causes statistical path length changes in the interferometer. In general, the phase change can be approximated as ΔΦ=*β*+*η(k-k*
_{0}
*)* with *k*
_{0} being the central wavenumber. It comprises a constant phase term *β* and a group dispersion term *η*=d(ΔΦ)/d*k* at *k*
_{0}. For the case of a single interface in the sample arm, moving with a constant axial velocity *V _{s}*, the phase shift can be written as ΔΦ=2

*kV*=2

_{S}τ*k*

_{0}Δ

*z*+2(

*k*-

*k*

_{0})Δ

*z*, i.e.

*β*=2

*k*

_{0}Δ

*z*and

*η*=2Δ

*z*. The change in optical path length Δ

*z*can also be introduced on purpose as in phase shifting FDOCT, using piezo actuators or electro-optic modulators (EOM). The latter case will involve higher order terms in

*k*for the resulting total phase shift due to the dispersion of the EOM crystal, which usually is Lithium Niobate (LiNbO

_{3}). An example for purely achromatic phase shifting, i.e.

*η*=0, has recently been given by employing acousto-optic frequency shifters to realize heterodyne FDOCT detection [29].

Keeping the linear dependence of the phase shift on *k*, we continue in calculating the FDOCT signal that is obtained by Fourier transforming Eq. (2) with respect to *k*′=*k*-*k*
_{0}. Expressing the sinc-function in Eq. (2) as
$\mathrm{sin}\mathrm{c}\left(\frac{\beta +\eta \left(k-{k}_{0}\right)}{2\pi}\right)=\mathrm{sin}\mathrm{c}\left(\frac{\mathit{\eta k\prime}}{2\pi}\right)\otimes \delta \left(\mathit{k\prime}+\frac{\beta}{\eta}\right)$
, and applying the convolution theorem we obtain

where $\widehat{\bullet}$
denotes the Fourier transformed signal. We see that for *β*=0, i.e. a change in group dispersion, the original axial point spread function (PSF), which is determined by the temporal coherence function of the light source, is averaged and broadened, due to the convolution with a *rect* function. However, for an exact signal analysis we need to take into account the complex exponential term exp(*jβz*/*η*) which is direct result from our spectral bandwidth centered at *k*
_{0}.

Assuming a Gaussian source spectrum with a full-width-at-half-maximum (FWHM) value Δ*k*
_{FWHM} centered at *k*
_{0} and writing the convolution in Eq. (3) explicitly one finds

with
$\alpha =\frac{\mathrm{\Delta}{k}_{\mathrm{FWHM}}}{\left(4\sqrt{\mathrm{ln}\left(2\right)}\right)}$
being the inverse of the round trip coherence length. The effect of the *rect*-function is to change the integration borders to *z*±*η*/2. The result of this integral will exhibit two signal peaks, centered at *z*=2*z*
_{0} and *z*=-2*z*
_{0}, symmetric around the zero path delay due to the complex ambiguity of the FDOCT signal. Keeping only the peak at *z*=2*z*
_{0}, and using iterative partial integration together with the definition of the Hermite polynomial of order *m*, H* _{m}*=(-1)

*exp(*

^{m}*x*

^{2})

*d*(exp(-

^{m}*x*

^{2}))/

*dx*, gives

^{m}The different summation terms in Eq. (5) are weighted by *αη*/*β*. In case of a moving sample interface we have
$\frac{\mathit{\alpha \eta}}{\mathit{\beta}}=\frac{\alpha}{{k}_{0}}=\frac{\mathrm{\Delta}{\lambda}_{\mathit{FWHM}}}{\left(2\sqrt{\mathit{ln}\left(2\right){\lambda}_{0}}\right)}$
, which is usually much smaller than one. In this case the contributions of *m*≥1 can be neglected. Upon using H_{0}=1, the absolute value of the expression in Eq. (5) can then be simplified to

Evaluating this expression at the interface position *z*
_{0}, i.e.(*z*-2*z*
_{0})=0 for a displacement during camera exposure Δ*z*, yields

Comparing this result with Eq. (2) shows the expected signal decay following a *sinc*-function. This is the result that would be expected for the achromatic case since it depends only on *β*. The effect of the *η*-term is an additional dampening factor that depends on the ratio between round trip coherence length 1/α and change in group dispersion *η*. For the moving sample interface, the first zero of the *sinc*-function is reached at Δ*z*=*λ*/2. Knowing the integration time *τ*, the corresponding sample velocity is then expressed as *V _{s}*=

*λ*/(2

*τ*). In the ideal achromatic case, the axial peak shape remains Gaussian. Plotting in Fig. 2 the axial peak shape for the moving sample interface as expressed by Eq. (6), it is revealed that the

*η*-term causes in addition peak broadening, which will eventually split up the signal and create two diverging signal peaks. The peak splitting is given by the width of the rect-function in Eq. (4) which in turn is equal to the displacement Δ

*z*. Nevertheless, as seen from Fig. (2), high SNR and large optical bandwidth is needed in order to observe this effect. Figure 3(a) indicates the normalized signal attenuation A(Δz), by plotting the normalized FDOCT signal peak heights [Eq. (6)] as a function of Δ

*z*.

#### 2.2. Signal recovery by reference phase tuning

In the preceding paragraph we mentioned the possibility to introduce phase shifting elements into the reference arm. If the phase shifting is performed during camera integration, then the overall phase shift ΔΦ [see Eq. (2)] will in general have contributions from both interferometer arms, i.e., ΔΦ=ΔΦ* _{R}*-ΔΦ

*. One immediately observes that for the phase matching condition*

_{S}the overall phase shift is zero, meaning that there is no interference fringe blurring and the signal will be completely recovered. The reference phase shifting can in principle be realized by a piezo actuator which is, however, relatively slow for fast and precise phase changes. We propose therefore to use an electro-optic modulator (EOM) in which case the phase change is realized by changing the refractive index of the EOM crystal fast axis. The corresponding phase shift will be a function of the voltage applied to the EOM. The relation *V*
_{π}(*k*)=(2*πa*-*bk*)/(*Qk*), the voltage at which for a certain wavenumber the phase is shifted by *π*, together with the constants *a*, *b* and *Q*, are provided by the manufacturer. The total change of the reference phase can be expressed as ΔΦ* _{R}*(

*k*)=

*πΔV*/

*V*, with Δ

_{π}*V*being the change of voltage during integration time. A Taylor expansion of

*V*

_{π}(

*k*) yields

*β*=

*πQk*

_{0}Δ

*V*/(2

*πa*-

*bk*

_{0}) and

*η*=2

*π*

^{2}

*Qa*Δ

*V*/(2

*πa*-

*bk*

_{0})

^{2}. It is then possible to associate a reference path length shift Δ

*z*

_{R}to the

*β*-term as

Finally, since the integration time is *τ*, we can write for the corresponding reference speed *V _{R}*=Δ

*z*/

_{R}*τ*, leading to a Doppler frequency of

*f*=

_{R}^{D}*k*

_{0}

*V*/

_{R}*π*. According to Eq. (8), this reference speed can compensate for the phase shift caused by sample movement, in particular of flow. For axial flow at constant speed

*V*we can write ΔΦ

_{S}*=2*

_{S}(k)*kV*. For

_{S}τ*V*=

_{S}*V*the

_{R}*β*-terms cancel out, but due to the group index of the EOM, the

*η*-terms will not compensate completely resulting in a signal peak broadening. Nevertheless, even for a large source bandwidth of 200

*nm*at a central wavelength of 800

*nm*, the difference to the ideally compensated situation will be smaller than 1

*dB*and can therefore be neglected. The important point is that the maximal detection sensitivity has been effectively shifted to the reference velocity

*V*, as shown in Fig. 3(b). It is worth noting that the maximum detectable velocity or Doppler frequency is independent of the Nyquist limit given by half the camera acquisition rate. The center Doppler frequency can in theory be shifted to arbitrary values only limited by the maximal EOM phase shift during camera exposure. It can be tuned to a particular velocity or associated Doppler frequency for which the sample Doppler frequency will be in resonance. The moving sample signal is then detected with high contrast whereas the signals of static structures are attenuated. The contrast between static and moving structures will be higher, the larger the actual sample velocity is.

_{R}#### 2.3. Differential velocity mapping

The freedom of driving the EOM with programmable phase slopes opens additional advantages of the resonant Doppler scheme. By driving the EOM with voltage ramps having inverse slopes, one obtains in effect two images where structures with opposite flow directions with respect to the optical axis appear enhanced. Flow inversion may be due to either different orientations of the same vessel or opposite flow direction in arteries and veins. A clear distinction is in fact possible within a true 3D reconstruction of the vessel bed. The two images can be subtracted having positive and negative intensity values in blue and red respectively to realize a bidirectional color Doppler flow map. Such representation has the additional advantage of removing spurious signals of static structures that will be visible depending on the reference velocity offset. One achieves a differential flow mapping already on an intensity basis without the need of evaluating FDOCT signal phase differences. It should be mentioned that the opposite reference velocities of course do not need to be the same; they can be individually tuned to the given sample flows in either direction by changing the respective EOM phase slope.

The above differential resonant Doppler scheme involves two phase slopes. We can also drive the EOM by one slope together with a constant voltage level. Knowing the exact reference velocity together with the theoretical signal decay, as shown in Fig. 4(a), one can then reconstruct the actual sample velocity based on the quotient between the two channels. Let us start with a reference velocity *V _{R}* and an a priori unknown sample velocity

*V*. Since we know the reference velocity we also know the normalized signal attenuation A(

_{S}*V*) and A(

_{S}*V*-

_{S}*V*). Having measured a signal level of

_{R}*I*̂

^{0}for constant EOM voltage and a signal

*I*̂

^{φ}for the channel shifted by

*V*, we then can write

_{R}and thus in principle

However, such reconstruction will be unambiguous only in the range [*V _{R}*-λ

_{0}/(2τ), λ

_{0}/(2τ)] indicated by the red shaded region in Fig. 4(a). Within this range one can associate a unique velocity value to a given signal attenuation difference ΔA as shown in Fig. 4(b). In practice this range is determined by the fact that a local signal can only be attenuated down to the noise level. This defines a maximum possible attenuation (MPA) (see Fig. 4). The velocity error is directly related to the SNR within the image channels

*I*̂

^{φ}and

*I*̂

^{0}. Writing for the logarithmic signal difference

*Q*(assuming equal noise

*δI*̂) with an error δ

*Q*, and defining

*I*̂

*as the theoretical signal level without blurring, we can write*

_{S}$$\phantom{\rule{3.em}{0ex}}\approx \mathit{log}\left({\mathit{\xce}}_{S}A\left({V}_{S}-{V}_{R}\right)+\mathit{\delta \xce}\right)-\mathit{log}\left({\mathit{\xce}}_{S}A\left({V}_{S}\right)+\mathit{\delta \xce}\right),$$

and we finally find for the statistical error of the differential image

where we used the definition *SNR*=(*I*̂* _{S}*/

*δI*̂)

^{2}and assumed that the signal intensities are large as compared to the noise. The expression in Eq. (13) can be compared to the phase difference error of phase-sensitive Doppler FDOCT being [12]

*δ(*Δ

*φ)*=1/√

*SNR*. The important difference is the dependence of the error on the signal attenuation for the static and the reference shifted channel A(

*V*) and A(

_{S}*V*-

_{S}*V*) respectively. The actual value of the absolute velocity error δ

_{R}*V*depends on δ

*Q*as well as on the local derivative of ΔA(

*V*) as shown in Fig. 4(b) as

_{S}The accessible unambiguous range and precision can be increased by using a three-stage EOM signal involving for example two inverse slopes together with a constant level, producing the tomograms *I*̂^{φ},*I*̂^{-φ} and *I*̂^{0} (see Fig. 5(A)). The only restriction is that the three attenuation curves need a common intersection part of their main lobes. A differential analysis of the three tomograms allows to determine the axial velocity of any sample motion such as flow in the range [-λ_{0}/(2τ), λ_{0}/(2τ)] (see §4). Note that the velocity bandwidth is dependent on the camera exposure time *τ* which is smaller than the camera line period T. The latter determines the unambiguous velocity bandwidth of phase-sensitive Doppler FDOCT as [-λ_{0}/(4τ), λ_{0}/(4τ)] [30̃].

An idea about the velocity can already be obtained by using the three different images as RGB channels for a color encoded Doppler image: red and blue encoding the opposite flow velocities and green displaying the static sample structure.

Finally it should be mentioned that, like in phase-sensitive Doppler FDOCT, only the axial component of flow *V _{ax}* is directly available. In principle a 3D reconstruction allows to extract the angle

*θ*of a vessel with respect to the optical axis. The actual speed would then be obtained from

*V*=

*V*/cos(

_{ax}*θ*).

## 3. Experimental

The setup is displayed in Fig. 5. It is a modified Mach-Zehnder configuration. The source is synthesized by two superluminescence diodes (EXALOS) with center wavelengths 827*nm* and 853*nm* and respective full-width-at-half-maximum values (FWHM) of 25*nm* and 34*nm* resulting in an overall FWHM of 36*nm*, a central wavelength of 833.5*nm* and an axial resolution of 8.5*μm* in air. The light is first linearly polarized (Pol) in order to match the fast axis orientation of the EOM’s LiNbO_{3} crystal. The reference arm contains the non-resonant EOM (NovaPhase) for phase shifting and a quarter waveplate in order to realize circular polarization. The EOM is driven by an arbitrary function generator (Agilent 33220A) followed by a high voltage amplifier (20x). The dispersion compensation matches both sample arm optics as well as the water chamber of the eye. The translation stage (TS) helps to adjust the reference to the sample arm length. The sample arm contains a LiNbO_{3} crystal to match the dispersion of the EOM in the reference arm. The light is then expanded with a Galilei telescope of focal lengths f_{1}=-40*mm* and f_{2}=100*mm*. The light passes the X-Y galvo scanning stage (X-Y Sc) (Cambridge Technology) and illuminates the eye via lenses f_{3}=60*mm* and f_{4}=30*mm*. The optical power at the cornea is 300*μW* which is safe for direct beam viewing according to the ANSI laser safety regulations [31]. After recombination, the reference and sample arm light is guided through a single mode fiber to the spectrometer module. The spectrometer is equipped with a volume diffraction grating (DG) (Wasatch, 1200lines/*mm*), and an objective lens (OL) with focal length 135*mm*. The CCD is a line scan camera (ATMEL AVIIVA, 2048 pixel, 12bit) where only 1024 elements are actually used. The camera is driven at a line rate of 17.4*kHz* with an integration time of *τ*=43*μs*. The full covered spectral width is 68*nm* resulting in a depth range in air of 2.6*mm*. The beam waist at the cornea is 1.8*mm*. The measured sensitivity is 98*dB* close to the zero delay with a sensitivity decay of -7*dB*/*mm*.

The synchronization of the EOM with the camera integration is shown in Fig. 5(A) showing the three-stage mode as explained at the end of §2.3. For a maximal voltage slope applied, the realized EOM phase shifts during the camera integration time of *τ*=43*μs* are ΔΦ* _{ref}(k)*=±2

*π*at the center wavelength of 833.5

*nm*corresponding to a velocity of

*V*=

_{R}*λ*/2

*τ*=±9.7

*mm*/

*s*and an associated Doppler frequency of ±23.3

*kHz*. The transverse scanning is performed continuously with a theoretical spot size of 14

*μm*at the retina.

## 4. Results and discussion

In a first step we measured the actual signal attenuation introduced by the EOM phase shifting by using a mirror as a sample. The curve showing the normalized signal decay with applied EOM voltage is plotted in Fig. 6 (circles). The measurement is in good agreement with the theoretical estimation according to Eq. (6) (solid line in Fig. 6).

This measurement is then repeated for an in-vivo measurement of static retinal structure close to the ONH. The movies in Fig. 7 show a time sequence at the same vertical position after correction of movement artifacts between tomograms using image registration [32, 33]. The EOM voltage slope was stepwise changed from 0.91 to 8.2*V*/*μs* in steps of 0.91*V*/*μs* during the recording of 81 tomograms. A single tomogram is acquired in 130*ms* and consists of 2250 lateral points covering a range of 3*mm* on the retina. Assuming the theoretic spot size of 14*μm* we thus have 10.5*x* over-sampling transversally. The recorded tomogram data has effectively three individual channels with 750 transverse points, each obtained quasi in parallel via the scheme outlined in Fig. 5(A): firstly, the standard FDOCT tomogram with attenuated flow signals due to interference fringe blurring [see Fig. 7(a)]; secondly, two tomograms where the static structure is gradually suppressed due to increasing EOM voltage slopes and the respective opposite flow signals are gradually enhanced [see Figs. 7(b) and 7(c)]. We chose a region in Figs. 7(b) and 7(c) where only static structure was present and calculated the average signal relative to the same region in the static channel [Fig. 7(a)] for each EOM voltage slope (triangles in Fig. 6). Each point represents an average over eight frames at the same voltage slope. As one can see, the measured structure attenuation deviates from the theoretical decay for higher voltage slopes. In order to explain this effect we calculated the MPA by referencing the average signal of the static structure within the indicated ROI in Fig. 7(a) to the average noise floor in the phase-shifted channels for each voltage slope (crosses in Fig. 6). The ROI was selected such that it contains only static structure avoiding flow signals. Along with the attenuation of the static structure the ratio of number of noise components to signal components within the ROIs of the phase-shifted channels increases. Thus, the signal attenuation will eventually reach the MPA value.

In the phase-shifted channels [Figs. 7(b), and 7(c)] one can clearly observe the loss of contrast for the static structures, whereas flow signals gain in contrast. The tomogram in Fig. 7(d) corresponds to the logarithmic differential signal between the two inversely phase-shifted channels, i.e. *log*(*I*̂^{φ}/*I*̂^{-φ}). In the logarithmic representation, the static structure is well suppressed whereas flow in opposite directions will have different signs and can easily be distinguished via appropriate adjustment of a grey level scale. This representation allows actually identifying roughly the velocity range by observing at which EOM voltage the flow signals will appear brightest or darkest. This in fact demonstrates clearly the idea of the resonant Doppler scheme. We observe that for vessels A and D [see Fig. 7(d)] the resonant region is situated at about ±6.4*V*/*μs*, corresponding to a flow velocity of approximately ±6.9*mm*/*s*, whereas for vessel B one finds a value at about -4.6*V*/*μs* with a corresponding velocity of -5*mm*/*s*. Only vessel C seems to gain gradually in contrast towards the maximum applied voltage slope of 8.2*V*/*μs*, hence *V*>8.9*mm*/*s*. The high flow velocities present in vessel C exceed in fact the detection bandwidth. Only within the bandwidth the signal will be visible. This bandwidth is now shifted with reference velocity during the time sequence. In addition we have flow changes because of blood pulsation during the data acquisition. Those effects together manifest as rings within the vessel that seem to move radially during the sequence. The velocity values are estimates of ±0.9*mm*/*s* given the voltage slope step size of 0.91*V*/*μs*. Comparing tomograms Figs. 7(b) and 7(d) to Fig. 7(a), the strong contrast enhancement realized by the introduced method is particularly well visible for vessels A and C. Since we do not correct for bulk motion velocity artifacts [12,17] they manifest as vertical band-like structures visible in Fig. 7(d).

The combination of the three channels as RGB entries of a color image results in the representation of Fig. 8. For low EOM voltage slopes the signal level will be approximately the same in all three channels resulting in an essentially black and white image. As the voltage increases the image turns gradually green, representing the static image channel, with the vessels tending towards blue and red (opposite flow directions). Hence the color within the vessels is also a qualitative indicator of the present velocities.

A clear quantitative value can be found by a differential analysis of all three channels as described in §2.3 and shown in Fig. 9. The tomograms are extracted from the time series in Fig. 7 at EOM voltage slopes of ±5.5*V*/*μs*. The left column shows the logarithmically scaled differential images *I*̂^{φ}/*I*̂^{-φ}, *I*̂^{φ}/*I*̂^{0} and *I*̂^{-φ}/*I*̂^{0} respectively with outlined vessels (A and B in Fig. 9). The right column shows the respective expected differential signal attenuation curves ΔA as function of flow velocity (cf. also Fig. 4). Combining the information of the three graphs, as is demonstrated for two selected vessels (A, B) in Fig. 9, one is able to find unique values for flow velocities via comparison and exclusion on a pure intensity basis. This can be easily realized by first observing from the upper row in Fig. 9 whether the corresponding velocity for a given intensity value is positive or negative. For a positive value one uses the red shaded range in the second row for unique velocity determination whereas for a negative value the blue shaded region in the third row. The maximum flow velocity which can be determined quantitatively is limited by the MPA value. We determined the MPA by calculating an average signal intensity and setting the MPA equal to the corresponding negative SNR, yielding MPA=-18.4*dB*. This in turn gives a range of quantifiable velocity of ±8.76*mm*/*s* (green dashed line in Fig. 9) which is not achievable with phase-sensitive Doppler FDOCT due to fringe blurring.

Figure 10 shows a corrected quantitative flow map of the tomogram from Fig. 9. It is now possible to extract flow profiles along the indicated cross sections through the four vessels A-D of Fig. 10. For speckle reduction an averaging filter with a 5*x*5-kernel was applied to the selected regions. The profile of vessel C has a flat top at 8.3*mm*/*s* since we hit the MPA limitation. The velocities in the central region of vessel C are too fast to be properly detected by at least two channels and quantitatively analyzed. The velocity calculation does not take into account any velocity offset due to proband motion. This explains the offset of the profile of vessel D with velocity peak value of -3.4*mm*/*s* and the sign change towards the edges of the vessel which is also visible for vessel A. The maximum velocity of vessel A is measured to be 5.2*mm*/*s*, whereas vessel D is at -5.1*mm*/*s*. According to Fig. 4(b), and assuming a typical SNR of 18.4*dB* in our tomogram with dynamic range of 39*dB*, the velocity error is ±550*μm*/*s*.

Applying vertical scanning, we finally recorded a full 3D data set of the ONH region where we have a typical loss of flow signals due to the presence of high flow velocities. Like before, the data set consists of 2250 transverse points covering an angle of 7° on the retina. Again we applied the three-stage EOM phase switching resulting in a threefold interlaced 3D data set. The voltage slopes were set to ±9.1*V*/*μs* in order to maximally contrast the vascular against the static structure for a 3D vessel sectioning. Each of the three sub-tomograms has 750 transverse points. A set of 81 tomograms was recorded vertically along a range of 4°. The 3D data set was recorded in 10.7*s*. As for the previous measurement, motion artifacts are corrected by image registration.

We calculated the z-projection of each 3D data channel and merged the resulting fundus images corresponding to the standard FDOCT image projection and the two reference phase shifted copies as RGB entries, resulting in Fig. 11(a). The RGB fundus image allows immediately distinction between opposite axial flow directions. A clear differentiation between arterial and venous blood flow can only be found by looking to the 3D reconstruction of the ONH vessel structure displayed in Fig. 11(b). The EOM allows not only dimming the static structure leaving only the vascular structure but also separating opposite flow. We reconstructed individual 3D images for both the positively and negatively shifted channels. The two 3D rendered structures where afterwards colored in red and blue respectively and merged into the 3D image sequence of Fig. 11(b). From this representation it is directly possible to identify points at which vessels change their orientation with respect to the optical axis.

Together with the color encoded directional information we have an unambiguous identification of an artery as the vessel that lead blood flow up from the nerve head channel to the retina, and the veins which transport the blood downwards the optical nerve (Fig. 11). For better contrast we combined both phase shifted channels to form a monochrome 3D image and displayed it in Fig. 12 as anaglyph stereo representation. This image, viewed through red-cyan filter goggles, gives a clear spatial impression of the vascular structure in the ONH region. Note also the fine vessel structure on the frontal side.

The method of intensity based flow detection and optical vessel vivisection works well for large sample velocities as is the case in the ONH region, since the static structure will be strongly attenuated due to the large shift of the signal attenuation curve A in Fig. 3(b). Nevertheless, this shift will also have an adverse effect on signals of flow at low speed. Tuning the EOM phase shift to these small flow velocities would result in only a weak static signal reduction, not sufficient anymore for an easy vessel sectioning. This is particularly the case if one approaches the foveal region.

## 5. Conclusion

In the present work we establish the theory of resonant Doppler imaging with spectrometer-based FDOCT and demonstrate its feasibility for the extraction of in-vivo retinal blood flow in 3D, i.e. optical vivisection of vascularization. The advantage of the resonant Doppler flow imaging modality is the possibility of achieving flow extraction without elaborate signal processing algorithms. Moreover we present a scheme of how to measure flow velocity values already on the intensity basis without the need to extract the signal phase. This method has an improved velocity detection range as compared to phase sensitive Doppler imaging. In addition the use of a programmable EOM has the flexibility, within its specifications, to tune the system sensitivity to arbitrary flow velocities independent of the CCD detector speed.

## Acknowledgments

We acknowledge the support of EXALOS AG (Switzerland) and financial support of the Swiss National Science Foundation (SNF grant No.205321-109704/1).

## References and links

**1. **R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express **11**,889–894 (2003). [CrossRef] [PubMed]

**2. **J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. **28**,2067–2069 (2003). [CrossRef] [PubMed]

**3. **N. Nassif, B. Cense, B. H. Park, S. H. Yun, T. C. Chen, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “In vivo human retinal imaging by ultrahigh-speed spectral domain optical coherence tomography,” Opt. Lett. **29**,480–482 (2004). [CrossRef] [PubMed]

**4. **G. Hausler and M. W. Lindner, “Coherence radar and spectral radar-new tools for dermatological diagnosis,” J. Biomed. Opt. **3**,21–31 (1998). [CrossRef]

**5. **A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. Elzaiat, “Measurement of Intraocular Distances by Backscattering Spectral Interferometry,” Opt. Comm. **117**,43–48 (1995). [CrossRef]

**6. **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**,457–463 (2002). [CrossRef] [PubMed]

**7. **R. Leitgeb, L. Schmetterer, W. Drexler, F. Berisha, C. Hitzenberger, M. Wojtkowski, T. Bajraszewski, and A. F. Fercher, “Real-time measurement of in-vitro and in-vivo blood flow with Fourier domain optical coherence tomography,” in *Coherence Domain Optical Methods And Optical Coherence Tomography In Biomedicine Viii*(2004), pp.141–146.

**8. **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**,1734–1746 (2005). [CrossRef] [PubMed]

**9. **U. Schmidt-Erfurth, R. A. Leitgeb, S. Michels, B. Povazay, S. Sacu, B. Hermann, C. Ahlers, H. Sattmann, C. Scholda, A. F. Fercher, and W. Drexler, “Three-dimensional ultrahigh-resolution optical coherence tomography of macular diseases,” Invest. Ophthalmol. Visual Sci. **46**,3393–3402 (2005). [CrossRef]

**10. **E. Gotzinger, M. Pircher, and C. K. Hitzenberger, “High speed spectral domain polarization sensitive optical coherence tomography of the human retina,” Opt. Express **13**,10217–10229 (2005). [CrossRef] [PubMed]

**11. **Y. Yasuno, S. Makita, Y. Sutoh, M. Itoh, and T. Yatagai, “Birefringence imaging of human skin by polarization-sensitive spectral interferometric optical coherence tomography,” Opt. Lett. **27**,1803–1805 (2002). [CrossRef]

**12. **B. H. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. J. Tearney, B. E. Bouma, and J. F. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 μm,” Opt. Express **13**,3931–3944 (2005). [CrossRef] [PubMed]

**13. **S. Makita, Y. Yasuno, T. Endo, M. Itoh, and T. Yatagai, “Polarization contrast imaging of biological tissues by polarization-sensitive Fourier-domain optical coherence tomography,” Appl. Opt. **45**,1142–1147 (2006). [CrossRef] [PubMed]

**14. **R. A. Leitgeb, L. Schmetterer, C. K. Hitzenberger, A. F. Fercher, F. Berisha, M. Wojtkowski, and T. Bajraszewski, “Real-time measurement of in vitro flow by Fourier-domain color Doppler optical coherence tomography,” Opt. Lett. **29**,171–173 (2004). [CrossRef] [PubMed]

**15. **R. A. Leitgeb, L. Schmetterer, W. Drexler, A. F. Fercher, R. J. Zawadzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express **11**,3116–3121 (2003). [CrossRef] [PubMed]

**16. **B. R. White, M. C. Pierce, N. Nassif, B. Cense, B. H. Park, G. J. Tearney, B. E. Bouma, T. C. Chen, and J. F. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical Doppler tomography,” Opt. Express **11**,3490–3497 (2003). [CrossRef] [PubMed]

**17. **S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express **14**,7821–7840 (2006). [CrossRef] [PubMed]

**18. **R. Leitgeb, M. Wojtkowski, A. Kowalczyk, C. K. Hitzenberger, M. Sticker, and A. F. Fercher, “Spectral measurement of absorption by spectroscopic frequency-domain optical coherence tomography,” Opt. Lett. **25**,820–822 (2000). [CrossRef]

**19. **R. A. Leitgeb, W. Drexler, B. Povazay, B. Hermann, H. Sattmann, and A. F. Fercher, “Spectroscopic Fourier Domain Optical Coherence Tomography: Principle, limitations, and applications,” in *Coherence Domain Optical Methods And Optical Coherence Tomography In Biomedicine Ix*(2005), pp.151–158.

**20. **J. Flammer, S. Orgul, V. P. Costa, N. Orzalesi, G. K. Krieglstein, L. M. Serra, V. X. Renard, and E. Stefansson, “The impact of ocular blood flow in glaucoma,” Prog. Retin. Eye Res. **21**,359–393 (2002). [CrossRef] [PubMed]

**21. **L. Schmetterer and M. Wolzt, “Ocular blood flow and associated functional deviations in diabetic retinopathy,” Diabetologia **42**,387–405 (1999). [CrossRef] [PubMed]

**22. **C. E. Riva, S. D. Cranstoun, J. E. Grunwald, and B. L. Petrig, “Choroidal blood flow in the foveal region of the human disc,” Invest. Ophthalmol. Visual Sci. **35**,4273–4281 (1994).

**23. **E. Friedman, “A hemodynamic model of the pathogenesis of age-related macular degeneration,” Am. J. Ophthalmol. **124**,677–682 (1997). [PubMed]

**24. **C. E. Riva, J. E. Grunwald, S. H. Sinclair, and B. L. Petrig, “Blood velocity and volumetric flow rate in human retinal vessels,” Invest. Ophthalmol. Visual Sci. **26**,1124–1132 (1985).

**25. **S. Yazdanfar, A. M. Rollins, and J. Izatt, “In vivo imaging of human retinal flow dynamics by color Doppleroptical coherence tomography,” Arch. Ophthalmol. **121**,235–239 (2003). [PubMed]

**26. **J. W. You, T. C. Chen, M. Mujat, B. H. Park, and J. F. de Boer, “Pulsed illumination spectral-domain optical coherence tomography for human retinal imaging,” Opt. Express **14**,6739–6748 (2006). [CrossRef] [PubMed]

**27. **J. K. Barton and S. Stromski, “Flow measurement without phase information in optical coherence tomography images,” Opt. Express **13**,5234–5239 (2005). [CrossRef] [PubMed]

**28. **S. H. Yun, G. J. Tearney, J. F. de Boer, and B. E. Bouma, “Motion artifacts in optical coherence tomography with frequency-domain ranging,” Opt. Express **12**,2977–2998 (2004). [CrossRef] [PubMed]

**29. **A. H. Bachmann, R. A. Leitgeb, and T. Lasser, “Heterodyne Fourier domain optical coherence tomography for full range probing with high axial resolution,” Opt. Express **14**,1487–1496 (2006). [CrossRef] [PubMed]

**30. **R. Leitgeb, L. Schmetterer, M. Wojtkowski, C. Hitzenberger, M. Sticker, and A. Fercher, “Flow Velocity Measurements by Frequency Domain Short Coherence Interferometry,” SPIE Proceedings **4619**,16–21 (2002). [CrossRef]

**31. **A. N. S. Institute, “American National Standards for Safe Use of Lasers, ANSI Z.136.1,” (2000).

**32. **P. Thevenaz, U. E. Ruttimann, and M. Unser, “A Pyramid Approach to Subpixel Registration Based on Intensity,” IEEE Transaction On Image Processing **7**,27–41 (1998). [CrossRef]

**33. **R. J. Zawadzki, S. M. Jones, S. S. Olivier, M. T. Zhao, B. A. Bower, J. A. Izatt, S. Choi, S. Laut, and J. S. Werner, “Adaptive-optics optical coherence tomography for high-resolution and high-speed 3D retinal in vivo imaging,” Opt. Express **13**,8532–8546 (2005). [CrossRef] [PubMed]