1.

Introduction

Optical coherence tomography (OCT) provides a few ways to measure relative motions and deformations. Speckle tracking¹ can reveal motions as small as the axial resolution of the OCT system typically near 5 μm. Smaller displacements can be measured using the interferometric nature of OCT. Phase differences between successive image frames (B-scans) detect axial motion of a few nanometers during the interval between acquisitions of successive frames.^2,3 With a few milliseconds between frames, phase measurements reveal axial velocities of a few micrometers per second. This phase-differencing technique has been applied to elastography of the eye ex vivo.^4,5 It has also been used in vivo to visualize the retinal vessels; in that application, those phase differences that remain consistent over the larger area of the image are considered bulk motion and removed numerically.⁶

In elastography, the interesting motions may vary smoothly across the imaged region. The B-scans are usually acquired with a transverse scan, so variations in bulk motion during the scan cause the appearance of a variation across the image. Therefore, bulk motion cannot be identified simply as the motion of the bulk portion of the image, and a different method of rejecting bulk motion is required.

A range of pulsatile motions are seen in the eye, from the spontaneous pulsation of veins seen when intraocular pressure exceeds venous diastolic pressure, to accommodative thickening of the lens by a few hundred micrometers per second,⁷ to the pulsatile swelling of the macula, only four micrometers in amplitude as measured by interferometry.⁸ In addition, under transient artificial changes in intraocular pressure, the optic disk flexes a few tens of micrometers.⁹

Partial-coherence interferometry can measure the pulsation of the retina of the human eye, which is only a few micrometers. The pulse amplitude was mapped as a function of lateral position on the retina,¹⁰ showing greater motion near the macula and especially near the optic disc, but without the guidance of an image, it is difficult to know exactly what is moving. These measurements used the cornea as the reference reflector, avoiding artifacts due to bulk motion of the eye.

Phase measurements are intrinsic in every spectral-domain OCT scan, but clinical imaging methods do not typically maintain a corneal reference reflection. If we can separate bulk motion from relative motions within the image data, we can measure elastography in vivo in the human eye.

2.

Data Acquisition

The OCT acquisition hardware has been described in Ref. 11. For some of the measurements reported here, we used an electro-optic modulator in the reference path, as described in Ref. 12, to provide phase shifts for full-range OCT reconstruction. The super-luminescent diode (Exalos, Inc., Schlieren, Switzerland) with 830-nm central wavelength has 45-nm full width at half maximum bandwidth, resulting in a depth resolution in air of 6.7 μm. The beam diameter at the cornea was 1.2 mm resulting in a nominal spot size of 22 μm on the retina. A complementary metal–oxide–semiconductor detector array (BASLER, Ahrensburg, Germany; Sprint spL4096-140 km) provided 100,000 axial scans (A-scans) per second.

The change in optical phase of the received light between successive frames tells us the axial displacement of the scattering tissue between acquisitions of those frames. It will be convenient to think in terms of the average axial velocity $v$ of the tissue between successive frames. The interesting differential motions in the retina are only a few micrometers per second. If the B-scan period is ${\tau }_{\mathrm{B}}$, the measured phase differences will be $2\times 1.36\times 2\pi v{\tau }_{\mathrm{B}}/\lambda$, where 1.36 approximates the group-velocity index of the tissue.

2.1.

Considerations of the Frame Rate

The detectable phase differences are limited by noise. Some retinal layers are only 10 dB above the noise floor, so the phase differences between measured returns from these layers have noise of $1/\sqrt{10}\approx 0.3$ rad [Eq. (2) of Ref. 13]. That noise, in terms of axial velocity of the tissue, is

Eq. (1)

$$\delta v\approx \frac{1}{\sqrt{10}}\frac{\lambda }{1.36·4\pi {\tau }_{\mathrm{B}}}=\frac{0.015\mu \text{m}}{{\tau }_{\mathrm{B}}}.$$

A longer time ${\tau }_{\mathrm{B}}$ between frames increases the measured phase differences and improves the signal-to-noise ratio. Eye motion, however, limits the useful frame period in two ways.

First, if the eye moves between frames a distance larger than the speckle size, successive B-scans lose their correlation. Correlated B-scans have phase differences determined by axial motion, but when successive B-scans lose their correlation, random phase gradients across the speckle pattern begin to contribute. (An analogous problem in Doppler OCT comes from lateral scanning, as explained in Sec. 3.3 of Ref. 13.) The data acquisition patterns used here repeat a single slice, so motion of the tissue perpendicular to the measurement plane completely changes the observed speckle pattern. The cooperative subjects imaged for the results presented here show stable speckle patterns for 100 ms and longer.

Second, time-varying bulk motion of the entire eye is more difficult to separate from the interesting spatial variation of motion within the frame if the frame rate is slow. When a patient is being imaged by a clinical OCT scanner, we typically see about $1\mathrm{mm}/\mathrm{s}$ axial motion of the posterior of the eye with a strong frequency component at the pulse rate, but also significant components up to 10 Hz. This system, like most clinical OCT systems, scans the imaging beam across the tissue to acquire each frame. Time-varying velocity of the eye as a whole during the transverse scan gives the appearance of a spatial variation in velocity.

Figure 1 shows the phase difference between successive frames that can arise from bulk motion of the eye

Eq. (2)

$$\frac{4\pi }{\lambda }\frac{1\mathrm{mm}}{\mathrm{s}}{\tau }_{\mathrm{B}}\sin (2\pi 10\mathrm{Hz}t),$$

for 10-Hz sinusoidal motion. When the sinusoidal amplitude is $1\mathrm{mm}/\mathrm{s}$, the eye moves several times the imaging wavelength within the time to acquire a frame, so there are several phase wraps. Figure 1 plots the expression above with these phase wraps. The vertical axis in Fig. 1 is scaled to express the measured phase in terms of axial velocity of tissue. The horizontal axis shows acquisition time $t$ marked with the times of acquisition of successive frames.

Fig. 1

Illustration of simplified bulk motion, sinusoidal motion of $1\mathrm{mm}/\mathrm{s}$ at 10 Hz, as it would appear as an additive component to the measured tissue velocity, across five frames with period ${\tau }_{\mathrm{B}}=5\mathrm{ms}$. With this frame period, the full range $\pm \pi$ of phase difference between frames corresponds to a range $\pm 30\mu \text{m}/\mathrm{s}$ in axial velocity.

The expected bulk motion is very large compared with the motion of interest and varies significantly across each B-scan, but it is continuous across B-scans. The interesting spatially varying motions are nearly periodic with the B-scans. If we limit our attention to the motions that are steady on the time scale of the frame rate, the interesting signal is found in a narrow frequency band near the B-scan rate $1/{\tau }_{\mathrm{B}}$ and its harmonics (Fig. 2). To the extent that $1/{\tau }_{\mathrm{B}}$ exceeds the bandwidth of interesting motion in the eye, we can filter out the bulk motion without affecting the interesting signal.

Fig. 2

Illustration of the spectral content of the axial motion in the sequence it is acquired. Solid lines: The power spectrum of the bulk motion observed in a tomogram sequence, and smooth curves indicating the expected frequency bands of spatially varying motion. Dashed line: Illustration of the low-pass filter used in Sec. 3. The frequency axis is marked at the B-scan rate $1/{\tau }_{\mathrm{B}}$ and its harmonics.

The bulk displacement of the eye has strong frequency components at the heart rate and its harmonics.¹⁴ From frame registration of a sequence of tomograms, we can measure the bulk motion of the eye and plot the power spectrum of that motion as the solid line in the center of Fig. 2. This particular power spectrum has strong peaks at 6 and 7 Hz. From a simultaneous recording of fingertip blood volume from a pulse oximeter, the heart rate during this data acquisition was 1.1 Hz.

A frame period ${\tau }_{\mathrm{B}}=5.0\mathrm{ms}$ sufficiently separates the frequency bands where we detect differential motions from the spectrum of bulk motion, reduces the noise estimate of Eq. (1) to $3\mu \text{m}/\mathrm{s}$, and was used for the results presented here. A 2.4-mm-wide tomogram was imaged repeatedly for 2 s.

3.

Data Processing

OCT images are composed of complex-valued pixels. We use $P(x,z,n)$ to denote the complex value corresponding to image coordinates $x$ and $z$ and frame number $n$. From successive pairs of frames, we form a complex cross-amplitude

The cross-amplitudes $C(x,z,n)$ typically show a strong phase ramp across each frame for the reasons illustrated in Fig. 1. Figure 3 shows the cross-amplitudes from a retinal image sequence with phase encoded as color. The sums of the cross-amplitudes along the axial direction provide a low-noise estimate of the phase due to the sum of local and bulk motions, also plotted in Fig. 3.

Fig. 3

An image of cross-amplitudes $C(x,z,n)$ with phase encoded as hue. The complex-valued sums $D(x,n)$ have a phase that corresponds to the axial velocities plotted below the image.

The phase of the column sums plotted in Fig. 3 shows significant drift from bulk motion, which we would like to low-pass filter and remove. The signal is often lost while scanning dark tissue, while the scanning mirrors are returning for the next scan, or when successive OCT frames are uncorrelated due to a saccade. Filtering the phase itself tends to introduce phase-wrap errors during these periods of loss of signal, so we filter the complex numbers holding that phase. The complex-valued axial sums of cross-amplitudes

The sequence of sums $D$ can be thought as a function of time, with $D(t)$ being the sum of the cross-amplitudes between A-scans acquired at times $t$ and $t+{\tau }_{\mathrm{B}}$. The complex products of successive sums,

The assumption implicit in applying a low-pass filter—that $S(t)$ is continuous across frame boundaries—is equivalent to the assumption of continuous bulk velocity between successive B-scans. Drop-outs in the data stream result in regions where the magnitude of $S(t)$ is small; the low-pass filter fills in these regions with values having phase near that on either side of the drop-outs. The filtering effectively interpolates the estimated bulk velocity across data drop-outs.

We implement the low-pass filter in the Fourier domain by multiplying with a raised cosine centered at zero frequency and reaching zero amplitude at the B-scan rate, as illustrated by the dashed curve in Fig. 2. The resulting estimate of bulk velocity, in terms of the resulting phase change $\mathrm{d}\mathrm{\Phi }$ between A-scans, is

There are discontinuous steps in phase at the boundaries between B-scans with seemingly random magnitudes on the order of 0.5 rad. These steps are likely due to the mechanics of the scanning mirrors, as they return to slightly different positions after each fly back. These steps in phase can be accounted for by rotating the phase $\mathrm{\Phi }$ in each frame to best match the average phase from the cross-amplitudes, $D(x,n)$, in that frame.

The estimated phase due to bulk motion is then removed from the cross-amplitude data

In order to improve the signal-to-noise, several frames are averaged together. The intensity frames $I$ are first registered using the StackReg¹⁵ plug-in for ImageJ.¹⁶ Then, the same registration shifts are applied to the corresponding velocity images $V$.

Even after image registration, the intensity $I$ at a given pixel can vary frame-to-frame. The cross-amplitude $C(x,z,n)$ can go to zero with corresponding random phase, if a saccade destroyed the correlation between the corresponding pair of frames. In order that the phases are properly weighted by their cross-correlation amplitudes, we average the cross-amplitudes

4.

Results

Acquisition on a quiescent section of the retina, an 8-deg wide scan placed a few degrees off the fovea, shows the degree to which bulk motion can be rejected. Figure 4 shows an average of 20 tomograms from a period free of saccades. The phase is quantified with contour lines spaced at $1\mu \text{m}/\mathrm{s}$, showing the average axial velocity of the tissue during acquisition of the 20 tomograms. Contours are drawn on spatially smoothed data; in this case, the cross-amplitudes $C(x,z,n)$ are smoothed with a 50-μm wide and 30-μm tall Gaussian kernel for the purpose of generating contours.

Fig. 4

Average of 20 cross-amplitude images overlaid with contour lines spaced at $1\mu \text{m}/\mathrm{s}$. Solid (green) contour lines indicate ascending tissue, dashed (red) contour lines indicate descending, and the zero-velocity contour is a thin black line.

The contours show a complex pattern of motion in the choroid. The apparent motion in the vitreous is spurious, because this portion of the image is at the noise floor and the phase noise is very large there. There is not more than $1\mu \text{m}/\mathrm{s}$ velocity within the neural retina.

Near the larger blood vessels, there are large relative motions correlated with heartbeat.¹⁷ Figure 5 shows results from two sessions imaging the same optic disc. Contour lines, every $3\mu \text{m}/\mathrm{s}$, quantify the axial motion. Motion around the large blood vessel, to the left of the cup, dominates.

Fig. 5

Averages of 20 tomograms from the optic disc overlaid with contour lines spaced at $3\mu \text{m}/\mathrm{s}$. This is a left eye, oriented conventionally with nasal side to the left and temporal to the right. The two images are from two imaging sessions of the same eye. The letters on the lower image indicate the locations of the data used for Fig. 6.

Very near to some blood vessels at the right of the image, there has been a phase wrap as the local velocities exceeded the $\pm 30\mu \text{m}/\mathrm{s}$ encoded as a full cycle of phase between tomograms, with the 5-ms between tomograms used here. The phase-wrap artifacts, however, only affect their local portions of the image.

At any point $(x,z)$ in the tomogram, we can extract the local tissue velocity as a function of time. Time traces from four locations, marked by letters a–d in Fig. 5, are plotted in Fig. 6. The curves have been smoothed by averaging the cross-amplitudes $C(x,z,t)$ over about 15 μm spatially and 50 ms temporally.

Fig. 6

Curves (a)–(d) are axial velocity as a function of measurement time (curves a–c being offset by multiples of $15\mu \text{m}/\mathrm{s}$) at the four locations indicated by the bulls-eye marks in Fig. 5: (a) the nerve fiber layer in the rim, just nasal of the large blood vessel lumen, (b) the cup, (c) the lamina cribrosa, and (d) the retinal pigment epithelium (RPE) at the temporal margin of the optic disk. Curve (e) is a sinusoid synchronous with the pulse, representing the Fourier component that is extracted from these and similar time traces to produce the image in Fig. 7.

Trace (a) shows motion of the nerve fiber layer near a large blood vessel along the optic-nerve rim. Trace (b) shows the large motion of the nerves at the bottom of the cup, probably caused by pulsation of blood vessels exiting the optic nerve head. Trace (c) from a visible portion of the lamina cribrosa shows much less motion. The retinal pigment epithelium (RPE) at the edge of the optic disk (d) shows motion counter to that of the cup.

The dominant motion in Fig. 6 is largely at the pulse frequency. We can extract from each point $(x,z)$ the Fourier component of the motion synchronous with the pulse and use these components to visualize the motion across the tomogram. These Fourier components, being amplitudes of the sinusoidal portion of the motion, are shown quantitatively as a contour plot in Fig. 7. The averaging of the data over the entire imaging time reduces velocity noise below $1\mu \text{m}/\mathrm{s}$, which is the spacing of the contour lines in Fig. 7.

Fig. 7

Fourier component of the motion synchronous with the pulse, found by Fourier analysis of all tomograms in a 2-s data stream. The contour lines are spaced at $1\mu \text{m}/\mathrm{s}$.

Figure 7 shows the anterior and left (nasal) portions of optic disc moving several micrometers relative to the posterior and right (temporal) portions synchronously with the heartbeat. We can see some of the lamina cribrosa in the center at the bottom of the image. The motion of the lamina at the pulse frequency has amplitude $4\mu \text{m}/\mathrm{s}$, but we see the RPE, too, on the right of the image moving with the same amplitude and direction. The flexing of the lamina relative to the optic disc opening appears to be only 1 or 2 μm.

5.

Conclusion

The phase difference between tomograms in an image sequence from an ophthalmic OCT system is dominated by the effects of bulk motion of the entire eye, but the bulk and differential motions can be practically separated. The acquisition rate should be sufficiently fast that speckle is correlated on successive tomograms. Steady differential motion within the tomogram is encoded in the data stream as a repetitive signal at the frame rate, so a fast frame rate is also beneficial for the frequency separation it provides between the differential and bulk motions. It is beneficial to work directly with the complex numbers that comprise the tomograms, and their products and sums, rather than their phases, because this allows proper averaging of the data in a natural way.

Relative velocities as small as $1\mu \text{m}/\mathrm{s}$ can be imaged in vivo in the human eye using frame-to-frame phase-sensitive OCT, despite the presence of much larger bulk motion.