2.3 Noise and error sources

The effect of shot noise on differential phase measurements has been described in detail by others [

22,

23]. In shot noise-limited conditions, the minimum measurable phase shift in the presence of photon noise (a zero-mean, complex Gaussian random process) is given by SNR

^{-0.5} [

24]. In addition to shot noise, there are additional factors that cause phase noise when imaging scattering tissue

*in vivo* with actively controlled galvanometers. These include physiologic motion (respiration, heart rate, etc.) and galvanometer jitter. The effects of these additional sources of noise can be divided into axial and transverse components. A shift in axial position causes both a phase shift and decorrelation. Phase shifts, caused by axial position shifts on the order of λ/4, may be corrected by motion correction algorithms [

25] since every axial position experiences the same phase shift. Decorrelation in the axial direction (caused by axial position shifts on the order of the axial resolution), leads to a different phase shift at axial positions separated by more than the axial resolution. Similarly, a shift in transverse position causes decorrelation [

26], also leading to a different phase shift at axial positions separated by more than the axial resolution. However, small amounts of positional noise predominantly cause axial phase shifts because an axial position shift of λ/4 leads to a phase shift of π, whereas position shifts must be on the order of the resolution (typically greater than several λ in both the axial and transverse directions) in order to cause significant decorrelation. It should be noted that algorithms exist to correct for axial phase shifts, while decorrelation caused by transverse motion is typically not corrected.

As shown in , galvanometer motion can affect phase measurements if the direction of the OCT scan direction does not lie in the plane of the wavefront of the incident beam. The OCT scan direction is defined as the direction of the three-dimensional displacement vector of the Gaussian beam during beam scanning, and is shown as a red arrow in . In general, the OCT scan direction forms a non-zero angle θ relative to the wavefront of the incident beam. The effects of this misalignment are twofold – 1) a projection of galvanometer angular “jitter” onto the axial direction, and 2) a Doppler shift as the OCT beam is scanned across static tissue, leading to a Doppler frequency shift from static tissue. The transverse displacement of the beam (v

_{x}T) is related to the incremental change in optical angle caused by the galvanometer (d

) by v

_{x}T =Ld

where L= (f

_{1}/f

_{2})f

_{OL}. As shown in

, any displacement along the axis of the scan is related to a displacement along the z axis by dz =sin(θ)v

_{x}T. Any z displacement causes a phase shift given by dΦ=(4πn/λ)dz, where n is the refractive index and λ is the wavelength. Therefore, the standard deviation of galvanometer jitter (σ

_{}) can be related to the standard deviation of phase fluctuations caused by the corresponding z displacement (σ

_{Φ}): σ

_{Φ}=(4πnL/λ)sin(θ)σ

_{}. Therefore, by minimizing the angle θ through careful optical alignment it is possible to minimize axial phase shifts caused by galvanometer jitter. In our setup we accommodated two steering mirrors in the beam path (M1 and M2), as shown in . The phase shifts along the x and y axes were monitored in real time, and the mirrors were adjusted to minimize the scanning-induced phase shift as the beam was scanned across a scattering phantom along each axis sequentially. Because both galvanometer mirrors cannot be simultaneously imaged at the objective back focal plane using our microscope design, it was impossible to achieve θ=0° at all transverse positions. However, this procedure achieved acceptably low residual Doppler shift levels for quantitative flow measurements. For qualitative angiography, correction of phase shifts due to bulk axial motion was performed as described in Section 2.7.

2.4 Bandwidth of static tissue and high-pass filtering

In order to understand the frequency sensitivity required for capillary visualization, we characterize the bandwidths of static and dynamic components of scattering from stationary tissue. For this purpose, we view the spatial speckle pattern generated by tissue as a complex white Gaussian random process filtered by the transverse point spread function [

27].

As both temporal and spatial relationships are important, we characterize the speckle pattern temporally and spatially by autocorrelation functions (R_{x}(x) and R_{t}(t)) and power spectral densities (P_{x}(f_{x}) and P_{t}(f_{t})). The variables f_{x} and f_{t} correspond to spatial and temporal frequency, respectively. The temporal and spatial autocorrelation functions are given by the Fourier transforms of the temporal and spatial power spectral densities.

The spatial and temporal autocorrelation functions have magnitudes with even symmetry about the origin and phases with odd symmetry about the origin.

Because the speckle pattern is generated by filtering of white Gaussian noise, its autocorrelation and power spectral density are determined solely by the power in the white Gaussian noise process and parameters of the transverse point spread function h(x). In order to calculate h(x), the transverse point spread function, it is necessary to calculate the electric field profile. Assuming Gaussian beam propagation, the incident electric field is given by:

For convenience we have normalized E

_{i}(x,y,z) so that the norm squared integrated over any xy plane is equal to 1. This ensures that the power in the Gaussian beam (proportional to the norm squared of the electric field) is independent of w

_{i}, the width of the beam waist.

As shown in , a point scatterer located at a position (x=x_{0},y=0,z=z_{0}) will be excited by the incident field and emit a spherical wave which is eventually coupled back to the single mode optical fiber. Therefore, the amplitude and phase of the emitted light is determined by the amplitude and phase of the incident electric field. In order to calculate the field amplitude and phase coupled back to the single mode fiber it is necessary to perform an overlap integral of the fiber mode and the backscattered field. While it is possible to perform this integration in any plane, it is most convenient to perform the integration in the plane of the point scatterer. The field of the point scatter in the plane z=z_{0} can be represented by

K is an arbitrary complex constant. The complex conjugate is necessary because the direction of propagation is reversed for the backscattered field. Therefore,

For convenience we assume that the point scatterer is located in the focal plane (z

_{0} = 0) and define the point spread function as:

Therefore, the point spread function is essentially given by the incident electric field squared.

Then, the autocorrelation function and power spectral density of the speckle pattern may be expressed as follows:

In the above expression σ

^{2} represents the variance (power) of white noise Gaussian random process and * denotes convolution. The power (variance) of the speckle pattern, R

_{x}(0), varies inversely as the beam width cubed. The fact that the OCT beam is scanned (x=v

_{x}t) imposes a relationship between spatial and temporal autocorrelation functions and power spectral densities.

The interpretation of this equation is that the temporal power spectral density of stationary tissue can be either stretched or compressed by increasing or decreasing the velocity, respectively. The variance of the temporal speckle pattern, or the integral of the power spectral density over all frequencies, is independent of v_{x}.

In order to verify these relationships we experimentally measured the spatial and temporal autocorrelation functions on a scattering phantom using different transverse scanning velocities. An OCT image of an Ocean Optics WS-1 diffuse reflectance standard is shown in
. A comparison of the measured and theoretical autocorrelation functions is shown in . The theoretical autocorrelation function is obtained from the measured Gaussian beam width and the theory presented above. shows measured spatial autocorrelation functions at different scanning speeds, demonstrating good agreement.

shows temporal autocorrelation functions at different scanning speeds, demonstrating an approximately twofold decrease in width as the speed is doubled.

, shows the spatial and temporal power spectral densities at different scanning speeds. The power spectral densities were obtained by Gaussian fitting of the autocorrelation functions followed by Fourier transformation. All measurements are normalized to the peak value of the slowest velocity (v_{x} = 0.55 mm/s) curve to allow easy comparison.

While shows the autocorrelation and power spectral density for stationary tissue, similar results hold for tissue moving in the axial direction. For tissue generating a Doppler frequency of f_{0}, assuming that the dominant source of decorrelation is transverse scanning, the magnitude of the autocorrelation remains unchanged while the phase is given by 2πf_{0}t. This leads to a shift of the power spectral density such that the maximum value is located at f_{0}. However, it is important to note that for moving red blood cells there will be additional sources of decorrelation leading to increased power spectral density width including axial speckle decorrelation caused by RBC motion in the axial direction, transverse speckle decorrelation caused by RBC motion in the transverse direction, and a random RBC velocity distribution.

In the rat somatosensory cortex, capillaries exhibit significant fluctuations in speed and flux, with speed varying between 0 and 1 mm / s and red blood cell flux varying between 0 and 100 s

^{−1} [

5]. Given the illumination and detection geometry, only velocities with components along the z axis generate Doppler shifts. Therefore, in order to map capillaries as well as large blood vessels, it is necessary to measure velocities much lower than 1 mm / s. For instance, a capillary velocity of 0.5 mm / s at an angle of 85 degrees to the probe beam generates a Doppler shift of approximately 150 Hz. In addition to Doppler frequency shifts caused by axial motion, broadening of the power spectral density is caused by speckle decorrelation. The effect of red blood cell motion (e.g. in capillaries) is to cause a speckle decorrelation, which leads to a broadening of the power spectral density. As a simple order-of-magnitude calculation of this effect, assuming a transverse motion of ~0.5 mm /s and referring to , the decorrelation time is on the order of 20 ms [], and the bandwidth generated is on the order of 50 Hz []. Broadening due to decorrelation from motion in the axial direction is also possible. The effect of power spectral density broadening is distinct from Doppler shifts, since it results in an increase in the width of the power spectral density as opposed to a shift in the center frequency of the power spectral density. In general red blood cell flow is accompanied by both frequency shifts and broadening.

The above discussion established estimates for the bandwidths of dynamic and static tissue under different conditions. The requirements for flow imaging and angiography are different; consequently they require different scanning protocols. Firstly, for flow quantification it is not necessary to measure capillary flow, since virtually all flow in the cortex can be traced to a surface artery or vein. In addition, repeated flow imaging may be advantageous to reduce physiologic flow noise. Therefore, a transverse scanning speed of ~21 mm / s was chosen for flow measurements, leading to a 3D OCT protocol of 160 images with 840 axial scans per image (protocol 1). This corresponds to an oversampling factor (defined as the number of A scans per FWHM of the autocorrelation function) of 10.5. Due to the short scanning time, this protocol was repeated 10 times in succession, reversing the direction of fast axis scanning for alternate 3D OCT scans. For angiography, it is desirable to reduce the temporal bandwidth of stationary tissue significantly to visualize capillary flows. A transverse scanning speed of 0.55 mm / s was chosen for angiography, leading to a 3D OCT protocol of 450 images with 32760 axial scans per image (protocol 2). This corresponds to an oversampling factor of 410. At this speed the bandwidth of stationary tissue was less than 50 Hz, as shown in .

A digital filter was designed to remove stationary scattering components from the OCT image [

16]. The digital filtering procedure assumes that the contribution from stationary and moving scatterers can be separated by a linear operation. (This is not the case if there is spatially correlated or multiple scattering, where the contribution of a stationary scatterer and a moving scatter in combination cannot be obtained by superposition of their individual contributions.) The ability of angiographic methods to visualize slow flow is determined by the bandwidths and magnitudes of stationary and dynamic scattering as well as the cutoff frequency of the high-pass filter. Because the amount of stationary scattering in the brain exceeded the amount of dynamic scattering, a relatively high cutoff frequency was chosen for the filters in
. A filter with a lower cutoff frequency will lead to a higher background from stationary tissue, but will enable imaging of slower flows, while a filter with a higher cutoff frequency will lead to a lower background from stationary tissue, but may also miss slow flows. An additional criterion in choosing the filter is that the filter should preserve the spatial resolution of the flow image. For protocol 1, the temporal extent of the filter was limited to ~0.2 milliseconds, corresponding to a spatial extent of ~4 µm. For protocol 2, the temporal extent of the filter was limited to ~3 milliseconds, corresponding to a spatial extent of ~1.7 µm and resulting is negligible loss of spatial resolution.

An additional consideration is the effect of high-pass filtering on power spectral density and frequency estimation.
shows the theoretical temporal power spectral densities before and after digital high-pass filtering for a range of Doppler shifts from 0 to 6 kHz. The width of the temporal power spectral density was assumed to be the same as that of static tissue for this figure. In other words, the only source of power spectral density broadening was assumed to be transverse scanning (Note that this assumption is not satisfied at high Doppler frequencies, corresponding to high axial velocities, due to additional decorrelation caused by axial motion). As can be seen by comparing power spectral densities before and after high-pass filtering in (protocol 1), high-pass filtering leads to two effects: a signal loss, and a distortion of the power spectral density. Because of this distortion, when the mean frequency is measured after high-pass filtering, a significant error results below 2 kHz []. By comparison, , shows that for protocol 2, there is negligible distortion and error in the mean frequency caused by high-pass filtering. Therefore, high-pass filtering leads to a distortion of the power spectral density, which can be mitigated by slower scanning speeds.

The preceding analysis does not account for continuous movement of the sample beam during the camera integration time. However, we note that the camera integration can be viewed as an analog low-pass filter, which can be analyzed in a similar manner to the high-pass filter in the preceding paragraph. This low-pass filter leads to two effects which have been described previously: a signal loss [

28] and an error in Doppler shift estimation [

29]. Both of these effects may be understood using the formalism we have described: At higher scanning speeds, the width of the temporal autocorrelation of static tissue is decreased, leading to increased temporal power spectral density width. If the temporal power spectral density width increases, static tissue experiences more attenuation from a temporal low-pass filter, leading to the signal loss described by Yun et. al [

28]. Additionally, if the temporal power spectral density is wide, low-pass filtering causes distortion and a frequency estimation error. This distortion should be worse at higher Doppler shift frequencies, as confirmed by Walther et. al [

29].

2.5 Velocity estimation

While angiograms can be used to visualize perfusion, quantitative interpretation of OCT angiograms is challenging. The brightness in an angiogram is related only to the amount of scattering from moving sources. Therefore, a decrease in intensity of the angiogram may be due to a decrease in scattering or a decrease in velocity, or a combination of both. Measurements of velocity and flow would therefore improve the utility of OCT for quantitative neuroscience research.

Doppler OCT measurements estimate either the phase slope of the autocorrelation function in the time domain or the center frequency of the power spectral density in the Fourier domain [

30]. Our procedure for velocity estimation is an extension of the Kasai autocorrelation algorithm [

31]. Velocity estimation was performed after high-pass filtering [

16]. The Kasai autocorrelation, shown below, has been used frequently in phase-resolved Doppler OCT measurements [

32]

In the above equation A

_{i} is the ith complex axial scan. This method can be viewed as a series of two steps. The first step is estimating the value of the autocorrelation function at a single time lag [

31].

The summation can be performed spatially in the axial and/or transverse direction over a window. Larger windows can be used to reduce noise at the expense of spatial resolution. The second step in the Kasai method is estimation of the phase derivative of the autocorrelation function at zero time lag, as shown below.

However, when sampling is dense, the autocorrelation function generally extends over more than two axial scans. Therefore, an extension of the standard approach may be desirable in the case of dense sampling. In this case, it is possible to calculate the complex autocorrelation function over a range of delays (T, 2T, etc.) and use phase information over a range of lags, as opposed to just a lag of one exposure time,

In this equation the autocorrelation is estimated at multiple delays. In our procedure, we first estimate the autocorrelation function as described and estimate the Doppler frequency f

_{t,est} as the frequency corresponding to the maximum value of the fft of the autocorrelation function estimate. This is equivalent to maximizing the power spectral density estimate:

A window of N = 7 was used, and a 256 point fft was used. The estimation was only performed where the estimate of R

_{t}(0) exceeded the noise threshold. The frequency estimate could then be related to the z projection of velocity using the following expression:

2.6 Flow calculation

The topography of blood supply to the cortex consists of pial arteries on the surface of the brain which branch to form diving arterioles supplying well-defined cortical capillary beds. The capillary beds are drained by ascending venules, with multiple venules draining a capillary bed supplied by a single arteriole. Ascending venules then converge to form pial veins. There are pathways for lateral and collateral flow in the cortical parenchyma, but, in general, flow in surface vessels should account for all flow in the parenchyma below [

33]. By measuring ascending venules at the locations where velocity measurements have the highest signal to noise ratio, it is possible to obtain precise measurements of blood flow. shows an

*en face* slice through the Doppler OCT volume at a cortical depth of approximately 50 µm. A pial artery (A) generates both positive and negative Doppler shifts as it undulates relative to the probe beam. Small arterioles (red) and venules (blue) are visible. Total blood flow measurements obtained by venous and arterial flow should theoretically be identical. Veins have less pulsatility, venules are more numerous than arterioles, and measured venous blood velocities fall within the Doppler detection range of our system. Therefore, we chose to measure total venous flow by measuring ascending venules directly before they joined pial vessels. Venules were identified based on the OCT angiogram as well as the

*en face* Doppler images.

A three-dimensional OCT data set consisting of 160 images with 840 axial scans per image over an 800 µm square area was acquired in approximately 7.5 seconds. The three-dimensional OCT data set, after Doppler processing, yields a three-dimensional map of the z-projection of velocity, v_{z}(x,y,z). We assume that the vector velocity field, **v**(x,y,z), has zero divergence and therefore represents a conservative field. Therefore, the surface integral of **v**(x,y,z) representing velocity flux, or flow in a vessel, is independent of the surface chosen, provided that the surface bisects the vessel.

Typically a surface in the xz or yz plane is chosen for this integral. This is the most obvious choice of plane as OCT images are typically acquired as cross-sections. For instance, in ophthalmology, a dual circumpapillary scanning protocol has been proposed, where two cross-sectional images are used to calculate the angle of the vessel relative to the optic axis [

34]. However, this choice of plane has the disadvantage that it requires explicit calculation of the angle of the vessel relative to the optic axis. shows a vessel inclined at an angle

with respect to the optic axis. Referring to this figure, it is possible to calculate the flow by integrating over the yz plane:

For this work, we obtain flow from a 3D OCT scan. We choose to perform the surface integration over the

*en face* (xy) plane, which is perpendicular to the measured velocity z component. shows the projection of this vessel in the

*en face* plane. shows an

*en face* cut through the 3D Doppler OCT volume. shows a zoom of the vessel marked by a white square. The area of this vessel in the xy plane is proportional to 1/|cos(

)|, as shown in , whereas the z projection of the velocity is proportional to cos(

). Therefore, these two effects cancel when the integration is performed in the

*en face* plane. Mathematically, by choosing the surface normal in the –z direction,

This procedure does not require explicit calculation of the vessel angle, since the surface normal is parallel to the component of velocity that is measured. Areas for integration were determined by manually choosing a depth and an area in the

*en face* plane for each vessel. The area for integration was chosen to be larger than the vessel cross-section.