## Abstract

Phase variance-based motion contrast is demonstrated using two phase analysis methods in a spectral domain optical coherence tomography system. Mobility contrast is demonstrated for an intensity matched Intralipid solution placed without flow within agarose wells. Vasculature oriented transversely to the imaging direction has been imaged for 3-4 dpf *in vivo* zebrafish using the phase variance contrast methods. 2D phase variance contrast images are demonstrated with imaging times only 25% higher than a Doppler flow image with comparable statistics. En face images created by integrating depth regions of 3D zebrafish intensity and phase variance contrast data demonstrate vasculature consistent with expected images.

©2007 Optical Society of America

## 1. Introduction

Spatial mapping of vascularization has important diagnostic benefits including but not limited to ocular diseases such as age-related macular degeneration [1]. Fluorescein angiography and indocyanine green angiography are currently the gold standards of vascular visualization techniques, with an injection supplying a fluorescent dye into blood vessels for imaging [2]. The invasiveness of these angiography techniques combined with possible complications [3,4] limits the screening capabilities for large populations. Optical coherence tomography (OCT) [5,6] is a non-invasive technique which can optically resolve three-dimensional structures with sensitivity to types of flow within the images.

Spectral domain optical coherence tomography (SDOCT)[7] is a high speed version of OCT which acquires an entire depth scan, called an A-scan, at the same time during the acquisition time of the spectrometer τ. For each A-scan, the intensity and phase of each depth resolved reflection can be calculated. With the phase information the most common form of motion contrast that can be calculated is the axial flow rate, defined as the flow component parallel to the imaging direction. The measurement of the axial flow rate v_{axial}(z), sometimes referred to as Doppler flow imaging, can be calculated using the average phase change at a given depth z for A-scans separated by time T:

In this definition, *λ* designates the mean wavelength of the illumination light source and n is the refractive index of the material in which the phase measurement was taken. Previously demonstrated flow imaging techniques with SDOCT use successive A-scans to measure phase changes such that *T* = *τ* in Eq. (1) [8–10]. The maximum axial flow component that can be quantitatively determined is limited by phase wrapping, which limits the maximum velocity component measurement to *v*
_{axial,max}(*z*) = ±*λ*/4*nτ*.

The Doppler flow contrast method only determines the axial flow component parallel to the imaging direction of the system. This component is designated by v(cosθ), where v is the velocity of the flow and θ is the angle of the flow direction relative to the imaging direction. For cases such as in the retina, the primary direction of the blood flow is nearly perpendicular to the imaging direction such that θ∼ 90°. In order to visualize the blood vessels in these cases with the Doppler technique, the velocity v of the flow must be large enough to have a sufficient axial component to surpass the phase noise at that point, which is determined by the local signal-to-noise ratio [11–13].

Phase variance is another form of phase contrast available with the information calculated with the SDOCT system. Methods demonstrating the use of phase variance in SDOCT to image vasculature have utilized the variance of phase changes over successive A-scans during transverse scanning [11,14]. The primary source of contrast from this type of method is the difference in phase noise created from transverse scanning across uncorrelated scatterers versus correlated scattering layers. Due to the structural nature of the phase variance measurement in this method an increase in phase contrast requires a sacrifice in transverse resolution of the contrast image, resulting in a limit to the smallest resolvable vascular structures.

We present two phase variance analysis techniques capable of identifying regions of contrast based on the motion of the scatterers without any degradation to the transverse resolution. These methods are demonstrated to visualize contrast between regions of static and mobile scatterers. Visualization of *in vivo* blood flow is also shown for vessel orientations perpendicular to the imaging direction.

## 2. Methods

Imaging demonstrated in this paper is performed with a SDOCT system, using a spectrometer capable of imaging with an A-scan rate of 25kHz with an integration time of 36.1μs. At the fastest acquisition rate, the measured system sensitivity is 97dB. Using a 840nm superluminescent diode with a 50nm bandwidth, the measured axial resolution of the system is 8.2μm in air, which corresponds to approximately 5.9μm in tissue. The transverse resolution measured at the illumination beam waist is approximately 4.4μm.

#### 2.1 Phase change limitations

The phase change measured for a given depth z with a time separation T, designated as Δϕ(z,T), has a combination of several factors affecting the accuracy of the measurement:

The motion of interest is from the scatterers located at the depth z, referred to as the phase change term Δϕ_{scatterer}(z,T). The total phase change Δϕ(z,T) also contains the bulk motion of the entire sample along the imaging (axial) direction Δϕ_{bulk}(T), caused by relative phase motion between the sample and the system and ideally is independent of depth z of the scatterers. Δϕ_{SNR}(z) designates a phase error associated with the local SNR of the data calculated at the depth *z* and is independent of time for a constant SNR. Experimental and theoretical results have determined that the standard deviation of SNR-limited phase error for phase changes has the form [11–13]:

Δϕ_{error,other}(z,T) encompasses the other phase errors which may occur for SDOCT phase measurements, including but not limited to phase changes caused by transverse motion across a scatterer at depth z [11], artifacts associated with limited depth sampling during axial motion of the sample or other computational based phase anomalies. By understanding the effects Δϕ_{bulk}(T), Δϕ_{SNR}(z), and Δϕ_{error,other}(z,T) have on the accuracy of phase measurements in a SDOCT system, improvements can be made to reduce the adverse effects of these terms in the phase contrast images.

#### 2.2 Bulk phase removal

During each imaging session, interferometer motion drift and sample motion creates a measurable bulk axial motion of the OCT image. Without any compensation, all phase measurements would be dominated by this relative motion between the system and sample. By using all of the image pixels to calculate the bulk motion contribution, pure noise data points can distort the calculations. In specific sample cases which contain a high reflection interface, this surface may be considered a stationary reference reflection to allow removal of the bulk axial motion. For most cases a reference reflector is not available or not stationary enough for phase removal of bulk motion. In these situations the phase changes measured from all of the sample depths can be used to calculate the axial motion, as long as the noise data effects can be reduced.

We used a weighted mean technique to calculate the bulk motion, utilizing phase changes from all depths within a chosen region of the A-scan. Each phase change calculated was weighted by the amplitude of the OCT signal for that depth to reduce the effect SNR-limited phase error σ_{Δϕ,SNR}(z) has on the calculation of the bulk motion removal [9]. For phase measurements separated by time T, the bulk axial motion is calculated as:

The signal amplitude I(z) used in Eq. (4) is related to the imaged OCT intensity, which is defined by 20 log(I(z)) in this system. The z summation for this equation is over a chosen depth region of the A-scan which contains the sample reflections. The phase changes Δϕ(z,T) used were conditioned to limit phase changes between -π and +π to avoid phase wrapping issues. With this method of bulk phase removal from the entire depth sample reflections, the accuracy of this method becomes limited to approximately the SNR-limited phase error of the strongest sample reflection in the entire depth. This places an expected lower limit on identifiable phase changes after compensating for the other effects.

To demonstrate the effect of bulk motion calculations without any intensity weighting method, phase changes from multiple depths within a paper sample are demonstrated in Fig. 1. With a bulk relative motion between the paper sample and system of approximately 1.5 radians, all of the phase changes over the depth of the sample were calculated. Figure 1(a) plots a histogram of all of the calculated phase changes, with no threshold applied regarding the OCT intensity for each data point. The phase changes calculated were originally limited to be between -π and +π to limit phase wrapping issues, but the 2π phase range was shifted by the calculated mean of the phase changes to reduce the effect of the low intensity terms on the statistical analysis. Figure 1(b) is a histogram of the same phase change data with a threshold applied to only display phase data from reflections with OCT intensity more than 15dB above the calculated noise level. By increasing the threshold level, the calculated standard deviation of the phase changes decreased as shown in Fig. 1(c). The upper threshold limit of the data was limited by the dynamic range of the sample reflections. Figure 1(c) also contains the maximum SNR-limited phase noise expected for the given threshold level of the data. It is noted that the calculated phase error for the sample motion deviates from the expected phase noise. The additional sources of error may be explained through computational artifacts, sample specific properties and transverse motion which occur during the measured axial motion. Further experimentation and analysis is required to reduce the effects of these additional noise sources.

#### 2.3 Phase variance contrast

With the compensation of the bulk sample motion, assuming no correlation between the various phase noise terms, the phase variance calculated from the corrected phase change Δϕ(z,T) - Δϕ_{bulk}(T) is determined by the summation of noise terms from Eq. (2):

Each of the variance terms correspond to the previously described phase changes except for σ^{2}
_{error,bulk}, which is the phase error created by the bulk motion subtraction method. The primary source of motion contrast information is σ^{2}
_{Δϕ,scatterer}(z,T), With many types of motion, the measured variance of motion increases with the time separation between phase measurements T. At least three types of motion can be observed using the phase variance contrast of the scatterers, including but not limited to: (i) the axial component of motion caused by Brownian-type random motion, (ii) phase effects due to transverse motion of the sample [11], and (iii) variations in the axial flow over the image acquisition time. The combination of all of these effects on the phase variance measurements will help to visualize locations of flow independent of the direction of flow relative to the imaging system. With the assumption of no time dependence on the main phase noise terms, increasing the time between phase measurements allows for an increase in the relative contribution of the motion variance signal σ^{2}
_{Δϕ,scatterer}(z,T) compared to the other phase noise components.

Two different scan modalities are presented to allow acquisition of 2D phase contrast images with a variety of time separations between phase measurements, which will be referred to as the MB-scan and BM-scan in this paper. Referencing OCT scanning terminology, an A-scan refers to a single reflectivity depth scan in the sample. A B-scan is the 2D cross sectional OCT image constructed from multiple A-scans taken over different transverse positions. An M-scan refers to multiple A-scans taken over time at the same transverse location. For the purpose of this paper, an MB-scan is a term used to describe multiple M-scans taken over different transverse positions to create a 2D OCT image of M-scans. A BM-scan on the other hand is multiple B-scans taken over time for the same transverse scan region. A graphical representation of a portion of a sample scan pattern for a MB-scan and a BM-scan is shown in Fig. 3.

#### 2.4 MB-scan

The purpose of the MB-scan is to acquire additional temporal information of the scatterer motion to improve motion contrast. Increasing the amount of A-scans taken within each M-scan not only increases the maximum time between phase images that can be measured, but the improved statistics can increase the dynamic range of standard Doppler flow imaging described in Eq. (1).

The basic contrast metric used for MB-scans with this system is *σ*
_{Δϕ}
^{2}(*z*,*T*
_{2})-*σ*
_{Δϕ}
^{2}(*z*,*T*
_{1}), so the phase noise terms described in Eq. (5) which are independent in time such as the SNR limited phase error can be cancelled. The remaining contrast terms expected are σ^{2}
_{Δϕ,scatterer}(z,T_{2})- σ^{2}
_{Δϕ,scatterer}(z,T_{1}), which usually requires T_{2}≫T_{1} for significant motion changes. For the example of Brownian motion of diffusing spheres, the expected form of the variance of scatterer motion is σ^{2}
_{Δϕ,scatterer}(z,T) = DT, where D is the diffusion constant. The calculated phase contrast *σ*
_{Δϕ}
^{2}(*z*,*T*
_{2})-*σ*
_{Δϕ}
^{2}(*z*,*T*
_{1}) in this case would be *D*(*T*
_{2}-*T*
_{1})≈*DT*
_{2} for T_{2}≫T_{1}.

Figure 4 demonstrates the visualization of scatterers undergoing Brownian motion without any flow introduced to the system. Figure 4(a) shows a schematic of a 2% agarose well which has been filled with 0.1% Intralipid, chosen in concentration to approximately match the agarose scattering OCT signal intensity. Figure 4(b) shows an averaged OCT intensity image of the filled agarose well as described in the schematic, demonstrating the limited intensity contrast. The highly reflective air/water interface in this sample appears very bright due to the 25dB dynamic range chosen for the image, saturating a large extent of the coherence function for this reflection. The phase contrast image using σ^{2}
_{Δϕ}(z,T_{2})-σ^{2}
_{Δϕ}(z,T_{1}) is plotted in Fig. 4(c) in the case of T_{2}=1ms and T_{1}=40μs, calculated from each time separation within the 200 A-scans acquired at each transverse location. The acquisition time of each M-scan collected is 8ms, resulting in a total data acquisition time of 1.6s for the image in this example. A threshold of *σ*
_{Δϕ}
^{2}(*z*,40*μs*)≤1 is applied to the contrast image, corresponding to an intensity threshold level SNR_{OCT}(*z*)≥1 as described by Eq. (3) if the phase variance between successive A-scans is dominated by SNR-limited phase noise.

Quantitative information about the mobility of the scatterers can be determined by the variance of the phase changes for a range of time separations. For the data used to create the images of Fig. 4, sub-linear diffusion is observed. It is possible that the concentration of Intralipid was high enough in this sample to cause crowding effects to alter the Brownian motion of the scatterers. To estimate the size of the scatterers, the phase variance calculated for the shortest time separations of the data was fit to the linear diffusion equation. The diffusional motion of the Intralipid in this case is roughly consistent with diffusing spheres in water with 200nm diameter. This is consistent with the expected diameter range of Intralipid scatterers.

The MB-scan is an extension of conventional OCT Doppler imaging techniques. Depending on the parameters chosen, this method can require significantly more time per transverse position compared to the original Doppler technique in order to get adequate motion contrast. The increased M-scan time limits the ultimate speed for 2D phase contrast imaging with this technique, which can not be improved through improvements to the A-scan acquisition speed. Observation of flow pulsatility is limited by the phase contrast image speed, which for demonstrated contrast images in this paper are 200 times slower than the acquisition speed of a single OCT B-scan of the same size.

#### 2.5 BM-scan

The BM-scan is intended to be an efficient acquisition method of phase contrast images by maximizing the time separation between phase measurements for a given 2D phase contrast image time. The expected motion contrast associated with flow and Brownian-type motion increases with an increased time separation between phase measurements. By using the phase differences between successive B-scans, the data throughput of phase contrast images is increased and the phase variance motion contrast can be increased as well. In order to ensure a constant time separation between phase measurements for each transverse location and maintain consistency of contrast across the image, successive B-scans are used for transverse scans of the same direction (every 2^{nd} image from bi-directional scanning or uni-directional scan with flyback).

The phase contrast σ^{2}
_{Δϕ}(z,T_{2})- σ^{2}
_{Δϕ}(z,T_{1}) used for MB-scans is not ideal for use in a BM-scan because the phase variance measured for phase changes between successive A-scans σ^{2}
_{Δϕ}(z,T_{1}) where T_{1}=τ contains additional phase errors (mentioned in Eq. (5) in the form of σ^{2}
_{error,other}(z)) which do not occur in the measurement of σ^{2}
_{Δϕ}(z,T_{2}). These additional errors are created due to the transverse scanning occurring between the phase measurements of successive A-scans [11]. Without using σ^{2}
_{Δϕ}(z,T_{1}) to remove the dominant SNR-limited noise term σ^{2}
_{Δϕ,SNR}(z), numerical estimates of the phase error can be made to compensate for this effect based on theory.

As derived in the appendix, the SNR-limited phase error is of the form of Eq (3):

The interferometric OCT signal for a reflection R(z) at a depth z is given by S(z)^{2} = (ητ)^{2}P_{R}P_{S}R(z)/(hν_{0})^{2}. The OCT noise N^{2} is described by the shot noise distribution where the mean <N^{2}> and the standard deviation σ_{N2} of the noise are equal and <N^{2}> =(ητ/hν_{0})P_{R}. The measured OCT signal Ĩ is a combination of the signal and noise components:

The signal amplitude *I* is the same as used in Eq. (4) for bulk motion calculation, with the OCT intensity defined as 20log(*I*) = 10log(|Ĩ|^{2}). For this theory, it is assumed at the noise amplitude N remains approximately constant relative the phase of the noise term *ϕ _{N}*, which is completely random. The magnitude of the OCT amplitude is calculated as:

Since the phase of the noise term is completely random, time averaging the magnitude of the measured intensity cancels cross term, which results in 〈|Ĩ|^{2}〉 = *S*
^{2} + 〈N^{2}〉. The OCT signal *S*
^{2} is calculated by removing the OCT noise such that *S*
^{2} = |〈|Ĩ|^{2}〉-〈N^{2}〉|. The absolute value function allows for log plotting of the signal data for noise terms below the expected mean noise level.

Estimation of the magnitude of the mean noise level <N^{2}> for the OCT system was achieved by fitting the phase variance data against OCT signal for the regime where S≫N and is negligibly affected by the bulk motion compensation phase error. For the data presented in Fig. 5, phase data for the OCT intensity 10 - 25 dB above the noise level was fit to Eq. (3). Figure 5 plots this estimated phase error along and the measured phase error from 200 successive A-scans as a function of the normalized OCT signal. Instead of using a mirror with a reflectance reduced by a neutral density filter as the sample, the data acquired used a paper sample which is comprised of a range of reflections over the entire depth.

For S^{2}/<N^{2}> ∼1, the discrepancy between the expected and measured phase error is most likely caused by the statistical variations of the noise, limiting the measurement of the interferometric signal. The measured phase error reaches a maximum level at approximately 3 radians^{2}, which is caused by the imposed limitations of maximum phase changes between -π and +π. For a completely random distribution of phase changes between -π and +π, the standard deviation is approximately 1.8 radians. This limits the expected maximum phase variance measured for a purely noise situation to approximately 3.2 radians^{2}. For the very high S^{2}/<N^{2}> case, the additional phase error introduced by the bulk motion removal algorithm has a non-negligible effect. Without using bulk motion removal, the measured phase error would be limited by the relative bulk motion between the system and sample.

## 3. Results and discussion

For the *in vivo* studies of vasculature imaging, 3-4 dpf (days post-fertilization) zebrafish (Danio rerio) were used. Zebrafish is ideal for imaging due to the fully developed vasculature system and the optical transparency of the animal. Figure 6 displays the anatomy and vasculature of a typical 3dpf zebrafish through confocal images and histological sections [15,16]. Figure 6(b) is an image produced of a zebrafish strain designed to express GFP (green fluorescent protein) in its vasculature for improved visualization. The lines appearing on Fig. 6(a) and (b) depict the approximate anatomical location where the histological section of (c) was taken from. The locations of the primary blood vessels of the zebrafish tail (dorsal aorta and axial vein) are labeled along with the major anatomical features on the histological image.

At time points earlier than 1 dpf, zebrafish embryos had phenylthiourea (PTU) applied to them, a chemical used routinely to inhibit pigment production in zebrafish. Without the pigment blocker, the embryo would increase the production of highly absorptive pigment regions which would interfere with imaging. After allowing the zebrafish to mature to the 3-4 dpf stage, the zebrafish to be imaged were anesthetized using a Tricaine solution and positioned within a low melt agarose solution to hold the animals stationary within agarose wells. The orientation of the zebrafish within these wells was on its side, allowing for confocal images similar to shown in Fig 6 to be produced. This alignment allowed for easier accessibility to transversely separate blood vessels in the tail, as well as orienting the flow in these vessels to be transverse to the imaging direction.

MB-scan data taken from 4dpf zebrafish imaged across the tail is presented in Fig. 7 at the same approximate location depicted by the scan lines in Fig. 6(a) and (b). All images shown were dimensionally scaled and cropped to highlight the region of interest. The presented MB-scan data is from 200 transverse locations over 1.6 s total acquisition time with acquisition over 100% of the scan cycle. There was no processing or filtering applied to the shown images beyond the previously described thresholds and data analysis. Figure 7(b) plots the phase variance contrast σ^{2}
_{Δϕ}(z,T_{2})- σ^{2}
_{Δϕ}(z,T_{1}) with T_{2}=1ms, T_{1}=40μs calculated from 200 A-scans at each transverse location, and the imaging parameters chosen are the same as for the Intralipid imaging of Fig. 4(c). The locations identified with the arrows in this image correspond to the expected vessels identified in Fig. 6(b) and (c), labeled as the dorsal aorta and the axial vein. Figures 7(c) and (d) plot the Doppler images of the average phase change over successive A-scans for 10 and 100 total A-scans, respectively with the same intensity threshold as in Fig. 7(b). The only expected difference between these images is an increased suppression of the SNR-limited phase error due to increased averaging. The scale of these flow images is +/- 200μm/s, which corresponds to a phase change of +/- 0.12 radians for time separations of τ=40μs. The phase scale for these images has been chosen to be approximately 4% of the maximum possible dynamic range of phase changes because it is an efficient way to demonstrate the improved vascular visualization of very low axial flow velocities of the Doppler method through increased averaging of phase change data. One of the main limitations of the Doppler images is the SNR-limited phase noise. For a Doppler image dynamic range of +/- 0.12 radians created using a total of 10 A-scans, Eq. (3) is used to determine that the expected phase noise will be within this range for SNR_{OCT}(*z*)≤9dB. Further reductions to the phase noise occur with increased averaging to the calculated phase changes. Thresholds have been placed on the phase variance and Doppler images at the OCT intensity such that the signal is equal to the mean noise level in the system.

The comparison of blood vessel contrast visualization is best demonstrated with Fig. 7(b) and (d). The Doppler image demonstrates the SNR-limiting phase noise suppression required to observe flow nearly perpendicular to the imaging direction. Figure 7(b) demonstrates the visualization capabilities of phase variance contrast, observing motion over a much larger spatial region than in the Doppler image.

The phase variance contrast σ^{2}
_{Δϕ}(z,T_{2})- σ^{2}
_{Δϕ}(z,T_{1}) imaged in Fig. 7(b) demonstrates phase change variance after the removal of SNR-limited phase noise. The value of this noise removal is demonstrated in Fig. 8, comparing phase variance σ^{2}
_{Δϕ}(z,T_{2}) with T_{2}=1ms to the previously presented contrast image of Fig. 7(b). A threshold is applied to both images corresponding to an OCT intensity level such that SNR_{OCT}(*z*) ≥ 1. The effect of the noise removal in this case is clearly observed by this comparison.

The BM-scan data used for Fig. 9 was taken over the same scan region on the zebrafish tail as the previously demonstrated MB-scan. In this case, with 200 transverse pixels and acquisition over 80% of the scan cycle, the phase contrast image uses a phase separation time T_{2}=10ms. Due to the increased time between phase measurements used in this contrast implementation, the maximum axial flow velocity that can be quantitatively measured with a Doppler method without phase wrapping is 250 times smaller than for successive A-scan phase changes in this system. To reduce the total imaging time required for the phase contrast image, 5 total B-scans were used to create the BM-scan data resulting in a total image acquisition time of 50ms. With the resulting reduced statistics for the phase contrast image, median filters were applied (3 pixels transverse, 3 pixels axially) to reduce the noise observed in the image.

Figure 9(b) plots the phase variance σ^{2}
_{Δϕ}(z,T) after the removal of the estimated phase noise described by Eq. (3), with a threshold applied based on the OCT intensity signal such that SNR_{OCT}(*z*)≤1. The arrows in Fig. 9(a) and (b) highlight areas corresponding to the dorsal aorta and axial vein regions identified using the MB-scan data in Fig. 7. These expected locations correspond to the contrast visualization occurring in the BM-scan phase contrast image. Additional phase variance contrast was also observed below the expected regions for the blood vessels. With the increased time between phase measurements as compared to the MB-scan, shadowing artifacts are now present below the regions of transverse blood flow due to the refractive index variations created within the vessels.

MB-scans and BM-scans have both demonstrated the ability to identify the motion associated with transverse flow not easily identified using Doppler OCT techniques. One important consideration in comparing techniques is the minimum acquisition speed required for each phase contrast image. For an MB-scan using an M-scan acquisition time of T_{0}, calculation of phase change variance requires at least four independent phase changes for adequate statistics, allowing for a maximum time separation of T_{0}/4 for phase changes. The demonstrated MB-scan technique created a 200 transverse pixel contrast image using M-scans composed of 200 A-scans for a total acquisition time of 1.6s, which results in a maximum time between phase changes of 2ms for proper statistics of phase variance. Reducing the time duration of each M-scan of the MB-scan would lower the total imaging time of the contrast image, but would also shorten the maximum time between phase changes that can be determined and reduce the motion contrast observed.

The BM-scan can decrease the total contrast imaging time by reducing the total number of B-scans used. By reducing the number of total B-scans to N, there are N-1 phase changes to be used for phase variance calculations. For the BM-scan method shown in Fig. 9, 4 phase changes were used to calculate the phase variance contrast. Limiting the total number of B-scans used will only limit the statistics used in variance calculations without reducing the motion contrast observed. Applying image filtering can be used to compensate for the lack of image smoothing effects due to decreased statistics. Faster 2D phase contrast image acquisition results in faster 3D phase contrast data acquisition, which is important to reduce the effects of sample motion.

Phase contrast 3D data sets were acquired using successive BM-scans across the heart of the zebrafish. Figure 10 presents the data through en face images created from 51 averaged OCT intensity and phase variance contrast scans taken over 480 μm x 280 μm. Figure 10(a) is created through the summation of the linear OCT intensity image over the depth of the zebrafish, presented in a logarithmic scale. The phase contrast en face image in Fig. 10(b) is created with a similar method, summing over the entire depth of the zebrafish and presented in a logarithmic scale. Figure 10(c) is a zoomed in version of the heart region of the GFP-labeled zebrafish confocal image presented in Fig. 6(b). The en face contrast image correlates with the expected vasculature confocal image for a similarly aged zebrafish.

Comparing the contrast capabilities of the MB-scan and the BM-scan, both are capable of visualizing transverse flow and mobility using phase variance measurements. While the MB-scan is capable of measuring a wide range of time separations of phase changes, the BM-scan is limited by the transverse scan speed capabilities of the imaging system. MB-scans are able to use Doppler imaging techniques with the same data that is able to monitor the progression of phase variance contrast over time. The increased flexibility of this method requires the sacrifice of 2D phase contrast imaging speed. BM-scans are very efficient in acquiring 2D phase contrast images which sacrifice flexibility and statistics for the benefit of imaging speed. Due to the reduction of the maximum measurable Doppler flow due to phase wrapping with a BM-scan, there are very few situations where the Doppler method would be applicable. One of the main limitations of the BM-scan is the susceptibility to sample motion during acquisition. With the increased time separation between B-scans is an increased possibility of measuring large displacements caused from sample motion. Additional work is required in this case to correct for the increased motion artifacts in the phase measurements. While MB-scans can measure SNR-limited phase noise for removal from the measured phase variance, BM-scans require numerical estimates to account for the same noise. Technological advances to improve A-scan acquisition rates would not improve 2D phase contrast imaging rates for MB-scans but would be able to increase BM-scan image sizes for the same phase contrast and imaging times. With limited data acquisition buffers, 3D contrast images are most efficiently created using multiple BM-scans.

## 4. Conclusion

We demonstrate two different phase variance contrast techniques capable of imaging regions of mobility with no flow as well as regions containing transverse flow. Phase variance contrast was presented from data sets images over zebrafish vasculature in the tail and the heart. MB-scans and BM-scans are able to image phase contrast with a variety of phase change time separations and speed of image acquisition. 3D phase contrast data sets are created using multiple BM-scans.

## Appendix

## Deriving SNR and phase noise for shot noise limited OCT performance

Define a spectral domain optical coherence tomography (SDOCT) system, where the power from reference arm arriving at the spectrometer is P_{R} and power from the sample arm arriving at the spectrometer is P_{S}. Assume that P_{S} ≪ P_{R}. Integration time of spectrometer is τ and the spectrometer contains M pixels used in k-space measurements. In terms of the number of electrons converted by the CCD of the spectrometer, the measurement in k-space is of the form:

*S _{j}*(

*k*) is the interferometric signal is defined as ${S}_{j}\left(k\right)=\eta \sqrt{{P}_{R}\left(k\right){P}_{S}\left(k\right){R}_{j}}\tau /hv$ for a sample reflection R

_{j}at optical path difference z

_{j}=2(z

_{Sj}- z

_{R}) of interferometer. The summation of this signal is taken over all of the sample reflections. Define ${P}_{R}=\sum _{k}{P}_{R}\left(k\right),\phantom{\rule{.2em}{0ex}}{P}_{S}=\sum _{k}{P}_{S}\left(k\right).$.

*F _{DC}*(

*k*) +

*N*(

*k*) combine to form the shot noise distribution of electrons. The mean number of electrons is given by

*F*(

_{DC}*k*) =

*ηP*(

_{R}*k*)

*τ*/

*hv*, where

*η*is the combined light collection and electron conversion efficiency of the spectrometer for photons of energy

*hv*.

*N*(

*k*) is the random portion of the Gaussian distribution with variance

*σ*

^{2}

_{N(k)}=

*ηP*(

_{R}*k*)

*τ*/

*hv*and zero mean. Using the property of Fourier transforms

*FT*(

*A*+

*B*) =

*FT*(

*A*) +

*FT*(

*B*), the Fourier transform of F(k) produces the OCT intensity amplitude

*Ĩ*(

*z*):

Breaking it into the real and imaginary components of the complex Fourier transform of M data points in k-space:

$$\phantom{\rule{2.6em}{0ex}}={I}_{\mathrm{Re}}\left(z\right)+i{I}_{\mathrm{Im}}\left(z\right)=\sum _{k}^{M}F\left(k\right)\mathrm{cos}\left(kz\right)-i\sum _{k}^{M}F\left(k\right)\mathrm{sin}\left(kz\right)$$

Assume the number of data points in FT, defined as M, is large enough such that a delta function accounts for the Fourier transform of an oscillatory function. Let the arbitrary choice of z_{j} ≥ 0 for all reflection locations be taken into account as well. These assumptions lead to:

## DC Term

$$=\delta \left(z\right)\frac{\eta {P}_{R}\tau}{hv}\phantom{\rule{12.5em}{0ex}}$$

## Interferometric Signal

Define the signal $\tilde{S}\left(z\right)=S\left(z\right)\mathrm{exp}\left(i{\varphi}_{S}\left(z\right)\right)=\sum _{k}^{M}\sum _{j}2{S}_{j}\left(k\right)\mathrm{cos}\left(k{z}_{j}+{\varphi}_{Sj}\right)\mathrm{exp}\left(-\mathrm{ikz}\right)$. Using the identity cos(*kz _{j}* +

*ϕ*) = cos(

_{Si}*kz*)cos(

_{j}*ϕ*)-sin(

_{Sj}*kz*)sin(

_{j}*ϕS*):

_{j}Assume *S _{j}*(

*k*) changes slowly compared to cos(

*kz*) and sin(

_{j}*kz*). With this assumption make the approximation for the summation: $\mathrm{cos}\left(kz\right)\mathrm{cos}\left(k{z}_{j}\right)=\mathrm{sin}\left(kz\right)\mathrm{sin}\left(k{z}_{j}\right)=\frac{1}{2}\delta \left(z-{z}_{j}\right)$

_{j}$$\phantom{\rule{4.8em}{0ex}}=\sum _{k}^{M}\sum _{j=1}{S}_{j}\left(k\right)\left(\mathrm{cos}\left({\varphi}_{Sj}\right)+i\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left({\varphi}_{Sj}\right)\right)\delta \left(z-{z}_{j}\right)=\sum _{j}\mathrm{exp}\left(i{\varphi}_{Sj}\right)\sum _{k}^{M}{S}_{j}\left(k\right)\delta \left(z-{z}_{j}\right)$$

$$=\frac{\eta \tau}{hv}\sum _{j}\mathrm{exp}\left(i{\varphi}_{Sj}\right)\sum _{k}^{M}\sqrt{{P}_{R}\left(k\right){P}_{S}\left(k\right){R}_{j}}\delta \left(z-{z}_{j}\right)\phantom{\rule{8em}{0ex}}$$

$$\phantom{\rule{.2em}{0ex}}=\frac{\eta \sqrt{{P}_{R}{P}_{S}\tau}}{hv}\sum _{j}\mathrm{exp}\left(i{\varphi}_{Sj}\right)\sqrt{{R}_{j}}\delta \left(z-{z}_{j}\right)=S\left(z\right)\mathrm{exp}\left(i{\varphi}_{S}\left(z\right)\right)\phantom{\rule{3em}{0ex}}$$

## Shot Noise Analysis

The noise calculated in OCT comes from the Fourier transform component of the noise distribution in k-space:

To understand how the noise transforms, Parseval’s theorem for finite length Fourier transforms is required.

Using $\tilde{f}\left(z\right)=\sum _{k}^{M}\phantom{\rule{.2em}{0ex}}f\left(k\right)\mathrm{exp}\left(-\mathrm{ikz}\right)$, where f(k) is a real function such that (*f*(*k*))* = *f*(*k*).

$$\phantom{\rule{13.6em}{0ex}}=\sum _{k}^{M}\phantom{\rule{.2em}{0ex}}f\left(k\right)\sum \phantom{\rule{.2em}{0ex}}_{z}^{M}\tilde{f}\left(z\right)\mathrm{exp}\left(\mathrm{ikz}\right)=\sum _{z}^{M}\sum \phantom{\rule{.2em}{0ex}}_{k}^{M}f\left(k\right)\mathrm{exp}\left(\mathrm{ikz}\right)\sum \phantom{\rule{.2em}{0ex}}_{k\prime}^{M}f\left(k\prime \right)\mathrm{exp}\left(-\mathrm{ik\prime z}\right)$$

$$\phantom{\rule{1.6em}{0ex}}=\sum _{k}^{M}\phantom{\rule{.2em}{0ex}}f\left(k\right)\sum _{k\prime}^{M}\phantom{\rule{.2em}{0ex}}f\left(k\prime \right)\sum _{z}^{M}\mathrm{exp}\left(i\left(k-k\prime \right)z\right)$$

$$\phantom{\rule{4.9em}{0ex}}=M\sum _{k}^{M}\phantom{\rule{.2em}{0ex}}f\left(k\right)\sum _{k\prime}^{M}f\left(k\prime \right)\delta \left(k-k\prime \right)=M\sum _{k}^{M}{\mid f\left(k\right)\mid}^{2}$$

This theorem is important for relating the Fourier transforms of the shot noise distribution in k-space N(k), measured on the CCD. With $\tilde{N}\left(z\right)=N\left(z\right)\mathrm{exp}\left(i{\varphi}_{N}\left(z\right)\right)=\sum _{k}^{M}N\left(k\right)\mathrm{exp}\left(\mathrm{ikz}\right)$

With the definition for the mean variation over all k-space measured as 〈*N*(*k*)〉 = 0 :

$${\u3008{\mid \tilde{N}\left(z\right)\mid}^{2}\u3009}_{z}=M{\u3008{{\sigma}^{2}}_{N\left(k\right)}\u3009}_{k}=\frac{\eta \tau}{hv}M{\u3008{P}_{R}\left(k\right)\u3009}_{k}=\frac{\eta {P}_{R}\tau}{hv}$$

Each component of Fourier Transform of random Gaussian noise distribution results in a Gaussian distribution as well. The real and imaginary components N_{Re}(z), N_{Im}(z) of the Fourier transform of the noise distribution N(k) are random Gaussian distributions, all centered around zero mean such that 〈*N*(*k*)〉 = 〈*N*
_{Re}(*z*)〉 = 〈*N*
_{Im}(*z*)〉 = 0. With each component being independent of each other, the phase of the noise *ϕ _{N}* (

*z*) is completely random. Determining the properties of the noise components:

$$\u3008{\mid \tilde{N}\left(z\right)\mid}^{2}\u3009=\u3008{N}_{\mathrm{Re}}{\left(z\right)}^{2}\u3009+\u3008{N}_{\mathrm{Im}}{\left(z\right)}^{2}\u3009={\sigma}_{N\mathrm{Re}}^{\phantom{\rule{1.2em}{0ex}}2}+{\sigma}_{N\phantom{\rule{.2em}{0ex}}\mathrm{Im}}^{\phantom{\rule{1.2em}{0ex}}2}$$

The transform components N_{Re}(z), N_{Im}(z) have identical distributions, which means that σ_{NRe}
^{2}=σ_{NIm}
^{2}. Therefore:

The probability distribution of the real noise component, which is identical to the imaginary component distribution, is calculated to be of the form:

The probability distribution of the noise amplitude N(z) is determined from the individual component distributions:

$$\phantom{\rule{.2em}{0ex}}=\underset{0}{\overset{N}{\int}}{{P}_{o}}^{2}\mathrm{exp}\left(-\left({N}^{2}\right)/2{\sigma}_{\mathrm{Nre}}^{\phantom{\rule{1.2em}{0ex}}2}\right)d{N}_{\mathrm{Re}}\phantom{\rule{0.9em}{0ex}}$$

Normalizing the distribution for the noise amplitude:

Using the probability distribution, the standard deviation of the magnitude of the noise N^{2} is calculated:

## OCT Calculations

Combining all of this analysis for z ≠ 0, the OCT intensity amplitude, given by the Fourier transform of F(k) is of the form:

$$\phantom{\rule{.2em}{0ex}}=S\left(z\right)\mathrm{exp}\left(i{\varphi}_{S}\left(z\right)\right)+N\left(z\right)\mathrm{exp}\left(i{\varphi}_{N}\left(z\right)\right)\phantom{\rule{1em}{0ex}}$$

The properties of the signal and noise terms have been derived previously. The magnitude of the OCT intensity can be calculated:

Averaging the OCT intensity, since *ϕ _{N}* is completely random:

## SNR Definition

System SNR sensitivity definition in OCT is described by the ratio of magnitude of the signal S^{2} where the reflection R=1 to the standard deviation of the noise magnitude:

## Definition of Phase Noise

With the OCT signal Ĩ(*z*), the calculated phase *ϕ*(*z*) can deviate from the expected sample phase *ϕ _{S}*(

*z*) depending on the relative noise properties. To determine the noise effects on the error on phase measurements, a probability analysis of the phase is required. Since the phase accuracy does not depend on the sample phase, set

*ϕ*(

_{S}*z*) = 0 for convenience. Using Eq. (27) in this case, the phase can be determined through trigonometry. The noise components N

_{Re}(z) and N

_{Im}(z) have the same Gaussian distribution described earlier in Eq. (23). For the case where S ≫ N, the phase determination can be simplified.

The probability distribution of the phase *ϕ*(*z*) in this case is calculated from the probability distribution of the noise using each OCT signal component S(z)^{2}:

which has a calculated variance of:

Phase error is dependant on the local signal to noise ratio for a given reflector. Phase changes measured for a reflector of OCT signal S^{2}(z) requires two phase measurements, each with phase error associated with it. The phase variance determined for phase changes is twice the value of the error for a single phase measurement.

## Acknowledgments

Authors would like to thank Carl Zeiss Meditec for their help in creating the imaging system. Further acknowledgements go out to Le Trinh and Julien Vermot for their help with the zebrafish work.

## References and links

**1. **S.O. Sykes, N.M. Bressler, M.G. Maguire, A.P. Schachat, and S.B. Bressler, “Detecting recurrent choroidal neovascularization. Comparison of clinical examination with and without fluorescein angiography,” Arch. Ophthalmology **112**, 1561 (1994). [CrossRef]

**2. **J.D. Gass, “Stereoscopic atlas of macular diseases,” 4th ed. (Mosby, 1997).

**3. **L.A. Yannuzzi, K.T. Rohrer, L.J. Tindel, R.S. Sobel, M.A. Costanza, W. Shields, and E. Zang, “Fluorescein angiography complication survey,” Ophthalmology **93**, 611 (1986). [PubMed]

**4. **M. Hope-Ross, L.A. Yannuzzi, E.S. Gragoudas, D.R. Guyer, J.S. Slakter, J.A. Sorenson, S. Krupsky, D.A. Orlock, and C.A. Puliafito, “Adverse reactions to indocyanine green,” Ophthalmology **101**, 529 (1994). [PubMed]

**5. **D. Huang, E.A. Swanson, C.P. Lin, J.S. Schuman, W.G. Stinson, W. Chang, M.R. Hee, T. Flotte, K. Gregory, C.A. Puliafito, and J.G Fujimoto, “Optical coherence tomography,” Science **254**, 1178 (1991). [CrossRef] [PubMed]

**6. **W. Drexler, U. Morgner, F.X. Kartner, C. Pitris, S.A. Boppart, X.D. Li, E.P. Ippen, and J.G Fujimoto, “In vivo ultrahigh resolution optical coherence tomography,” Opt. Lett. **24**, 1221 (1999). [CrossRef]

**7. **A.F. Fercher, C.K. Hitzenberger, G. Kamp, and S.Y Elzaiat, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. **117**, 43 (1995). [CrossRef]

**8. **A.F. Fercher, W. Drexler, C.K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. **66**, 239 (2003). [CrossRef]

**9. **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 coherence tomography,” Opt. Express **11**, 3490 (2003). [CrossRef] [PubMed]

**10. **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 (2004). [CrossRef] [PubMed]

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

**12. **S. Yazdanfar, C.H. Yang, M.V. Sarunic, and J.A. Izatt, “Frequency estimation precision in Doppler optical coherence tomography using the Cramer-Rao lower bound,” Opt. Express **13**, 410 (2005). [CrossRef] [PubMed]

**13. **B. Vakoc, S.H. Yun, J.F. de Boer, G.J. Tearney, and B.E. Bouma, “Phase-resolved optical frequency domain imaging,” Opt. Express **13**, 5483 (2005). [CrossRef] [PubMed]

**14. **S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical Coherence Angiography,” Opt. Express **14**, 7821 (2006). [CrossRef] [PubMed]

**15. **
“The Zebrafish Information Network”, www.zfin.org.

**16. **S. Isogai, M. Horiguchi, and B.M. Weinstein, “The vascular anatomy of the developing zebrafish: an atlas of embryonic and early larval development,” Developmental Biology **230**, 278 (2001). [CrossRef] [PubMed]