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.

 

Fig. 1. Typical high resolution retinal tomogram with “empty” blood vessels (^*) in the optic nerve head region. Tomogram size: 12mm x 1mm (lateral x depth) with 8.5μm axial resolution (in air).

Download Full Size | PDF

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_S(k) are expressed as functions of the wavenumber 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 ₀ 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_AC(k) and assuming a phase shift ϕ(t) that is linear in time, the integration yields [28]

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 ₀ ) with k ₀ being the central wavenumber. It comprises a constant phase term β and a group dispersion term η=d(ΔΦ)/dk at k ₀. 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 ΔΦ=2kV_Sτ=2k ₀Δz+2(k-k ₀)Δz, i.e. β=2k ₀Δ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₃). 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 ₀. Expressing the sinc-function in Eq. (2) as $\sin \mathrm{c}\left(\frac{\beta +\eta \left(k-k_0\right)}{2\pi }\right)=\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 ₀.

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

with $\alpha =\frac{\mathrm{\Delta }{k}_{\mathrm{FWHM}}}{\left(4\sqrt{\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=2z ₀ and z=-2z ₀, symmetric around the zero path delay due to the complex ambiguity of the FDOCT signal. Keeping only the peak at z=2z ₀, and using iterative partial integration together with the definition of the Hermite polynomial of order m, H_m=(-1)^mexp(x ²)d^m(exp(-x ²))/dx^m, gives

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₀=1, the absolute value of the expression in Eq. (5) can then be simplified to

Evaluating this expression at the interface position z ₀, i.e.(z-2z ₀)=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.

 

Fig. 2. Normalized signal attenuation A as a function of Δz, the displacement of a sample surface during the integration time, and z, the axial coordinate in dB. The three plots correspond to a Δk _FWHM of 0.2/μm, 1/μm and 2/μm. For larger spectral bandwidths, the signal decays faster with Δz and the splitting occurs earlier.

Download Full Size | PDF

 

Fig. 3. (a). Maximum value in z-direction of the normalized signal attenuation A as a function of displacement Δz for the same parameters as in Fig. 2. (b) A as a function of sample velocity for the static case (solid line) and with a reference velocity offset V_R (dashed line). As indicated, the reference velocity offset reduces static structure intensity (green arrow) and enhances moving structure intensity (red arrow).

Download Full Size | PDF

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-ΔΦ_S. One immediately observes that for the phase matching condition

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 ₀ΔV/(2πa-bk ₀) and η=2π ² QaΔV/(2πa-bk ₀)². 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 ₀ 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_S we can write ΔΦ_S(k)=2kV_Sτ. For V_S=V_R the β-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 200nm at a central wavelength of 800nm, the difference to the ideally compensated situation will be smaller than 1dB and can therefore be neglected. The important point is that the maximal detection sensitivity has been effectively shifted to the reference velocity V_R, 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.

 

Fig. 4. (a). Signal attenuation A for the static (solid line) and shifted case (dashed line). The red shaded region indicates the intersection of the main lobes. The green dash-dotted line indicates the MPA for a typical local signal in the tomogram. (b) Velocity as a function of the quotient ΔA of the shifted and the static case within the red shaded area (red solid line). The curve allows associating to a given ΔA value a unique sample velocity V_S. The velocity error δV [Eq. (14)] can be read from the dashed line. Within the green dash-dotted lines it is possible attributing a unique velocity value (see text for details).

Download Full Size | PDF

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_S. Since we know the reference velocity we also know the normalized signal attenuation A(V_S) and A(V_S-V_R). Having measured a signal level of Î⁰ for constant EOM voltage and a signal Î^φ for the channel shifted by V_R, we then can write

and thus in principle

However, such reconstruction will be unambiguous only in the range [V_R-λ₀/(2τ), λ₀/(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 Î^φ and Î⁰. Writing for the logarithmic signal difference Q (assuming equal noise δÎ) with an error δQ, and defining Î_S as the theoretical signal level without blurring, we can write

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

where we used the definition SNR=(Î_S/δÎ)² 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_S) and A(V_S-V_R) respectively. The actual value of the absolute velocity error δV depends on δQ as well as on the local derivative of ΔA(V_S) as shown in Fig. 4(b) as

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 Î^φ,Î^-φ and Î⁰ (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 [-λ₀/(2τ), λ₀/(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 [-λ₀/(4τ), λ₀/(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_ax/cos(θ).

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 827nm and 853nm and respective full-width-at-half-maximum values (FWHM) of 25nm and 34nm resulting in an overall FWHM of 36nm, a central wavelength of 833.5nm 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₃ 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₃ 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₁=-40mm and f₂=100mm. The light passes the X-Y galvo scanning stage (X-Y Sc) (Cambridge Technology) and illuminates the eye via lenses f₃=60mm and f₄=30mm. 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 135mm. 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.4kHz with an integration time of τ=43μs. The full covered spectral width is 68nm resulting in a depth range in air of 2.6mm. The beam waist at the cornea is 1.8mm. The measured sensitivity is 98dB close to the zero delay with a sensitivity decay of -7dB/mm.

 

Fig. 5. Optical setup (see text for details).

Download Full Size | PDF

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.5nm corresponding to a velocity of V_R=λ/2τ=±9.7mm/s and an associated Doppler frequency of ±23.3kHz. 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).

 

Fig. 6. Normalized signal attenuation A for different samples due to the applied EOM phase slope: mirror (circles), in-vivo retinal structure (triangles) for ROIs in Figs. 7(a)-7(c) compared to theory (solid line) [Eq. (6)]. Crosses represent the maximum possible attenuation (MPA) of the static structure in Fig. 7(a).

Download Full Size | PDF

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.2V/μs in steps of 0.91V/μs during the recording of 81 tomograms. A single tomogram is acquired in 130ms and consists of 2250 lateral points covering a range of 3mm on the retina. Assuming the theoretic spot size of 14μm we thus have 10.5x 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.

 

Fig. 7. Time sequence (after registration) with a three-stage EOM signal applied showing the same vertical position on the retina with stepwise increasing EOM voltage; all tomograms in log-scale. (a) (1.8 MB) Static retinal structure which does not change with increasing phase shift applied,[Media 1] (b,c) (1.7 MB, 1.7 MB) phase-shifted tomograms with opposite phase shifting directions and [Media 2][Media 3](d) (1.9 MB) calculated tomogram representing Î^φ/Î^-φ. Axial flow directions are clearly visible for all four vessels. Tomogram size: 2.75mm x 1.24mm (lateral x depth) with 8.5μm axial resolution (in air). [Media 4]

Download Full Size | PDF

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(Î^φ/Î^-φ). 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.4V/μs, corresponding to a flow velocity of approximately ±6.9mm/s, whereas for vessel B one finds a value at about -4.6V/μs with a corresponding velocity of -5mm/s. Only vessel C seems to gain gradually in contrast towards the maximum applied voltage slope of 8.2V/μs, hence V>8.9mm/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.9mm/s given the voltage slope step size of 0.91V/μ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).

 

Fig. 8. (1.8 MB) Time sequence of Fig. 7(a)-7(c) in RGB representation. [Media 5]

Download Full Size | PDF

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.5V/μs. The left column shows the logarithmically scaled differential images Î^φ/Î^-φ, Î^φ/Î⁰ and Î^-φ/Î⁰ 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.4dB. This in turn gives a range of quantifiable velocity of ±8.76mm/s (green dashed line in Fig. 9) which is not achievable with phase-sensitive Doppler FDOCT due to fringe blurring.

 

Fig. 9. Velocity determination by unambiguous differential velocity mapping (see §2.3). For each differential image (lhs) a corresponding theoretical ΔA curve (rhs) is shown. The blue and red solid lines correspond to the selected vessels on the lhs. The horizontal lines and circles on the rhs indicate possible velocity values according to the intensity level in the corresponding tomogram. The ambiguity is removed by looking for the common velocity in all three differential maps as indicated by the respective vertical lines on the rhs. The green dashed lines indicate the MPA limit of -18.4dB. Tomogram size: 3mm x 1.65mm (lateral x depth) with 8.5μm axial resolution (in air).

Download Full Size | PDF

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 5x5-kernel was applied to the selected regions. The profile of vessel C has a flat top at 8.3mm/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.4mm/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.2mm/s, whereas vessel D is at -5.1mm/s. According to Fig. 4(b), and assuming a typical SNR of 18.4dB in our tomogram with dynamic range of 39dB, 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.1V/μ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.7s. As for the previous measurement, motion artifacts are corrected by image registration.

 

Fig. 10. (lhs) Velocity map obtained by differential analysis of Fig. 9. with zoomed and average-filtered ROIs of vessel regions. (rhs) Velocity profiles extracted along the indicated lines in the vessel regions.

Download Full Size | PDF

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.

 

Fig. 11. (a). RGB fundus image showing different axial flow directions in blue and red, and the static structure in green. Size: 2.4mm x 1.7mm. (b) (2.5 MB) Movie of merged 3D volumes of positively and negatively shifted data sets in red and blue respectively. Tomogram size: 2.4mm x 1.7mm x 1.1mm (length x width x depth) with 8.5μm axial resolution (in air). [Media 6]

Download Full Size | PDF

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.

 

Fig. 12. (3.6 MB) ONH blood vessel structure of Fig. 11 in an anaglyph stereo representation (red-cyan goggles). [Media 7]

Download Full Size | PDF

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.