Optical coherence tomography (OCT) has been used to measure micro-displacements of a biological specimen. This displacement measurement has been used to characterize tissue biomechanical properties [1] and monitor the dynamic change of retinal tissue during laser coagulation [2,3].

Several methods have been proposed to calculate this displacement measurement. Speckle tracking is a commonly used technique that measures displacements by finding maximum cross-correlation coefficients [1]. Phase-sensitive measurement is another approach that provides nano-scale displacements by analyzing the phase of OCT signals [2,4,5]. Recently, the value of correlation coefficients has also been used to measure the micro-displacements of tissues [3,6].

Among these methods, Kurokawa et al. proposed a method to measure lateral and axial displacements by simultaneously using Doppler phase shifts and spatially localized correlation coefficients. The method successfully generated dynamic micro-displacement maps of retinal tissue during laser coagulation [3]. However, the method does not provide the direction of the lateral displacement. In addition, the approach is based on the assumption that the sample displaces only within the imaging plane of an OCT B-scan. Technically, this method cannot distinguish between an out-of-plane displacement and an in-plane displacement.

In this Letter, we propose a new method to obtain spatially resolved in-plane bidirectional displacement and out-of-plane displacement between two OCT B-scan images using localized correlation coefficients.

Assuming a probe beam with a Gaussian profile on the sample and a Gaussian temporal coherence function for the light source, the point spread function (PSF) of OCT can be expressed as $h(\mathbf{r})=\exp [-2(x^2+y^2)/w_l^2-2z^2/w_z^2]\exp [-2i{k}_c(z-z_0)]$, where $x$ is the in-plane lateral position, $y$ is the out-of-plane lateral position, and $z$ is the depth position. Further, $w_l$ is the lateral $1/e^2$ spot radius of the Gaussian probe beam, $w_z$ is the axial $1/e^2$ radius of the PSF, $z_0$ is the zero-delay depth position of the OCT interferometer, and $k_c$ is the center wavenumber of the light source. Here, we consider the local displacement of a sample at position $\mathbf{r}$; $\mathrm{\Delta }\mathbf{r}(\mathbf{r})=(\mathrm{\Delta }x(\mathbf{r}),\mathrm{\Delta }y(\mathbf{r}),\mathrm{\Delta }z(\mathbf{r}))$, where $\mathbf{r}=(x,y,z)$ and $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, and $\mathrm{\Delta }z$ are displacements along the $x$, $y$, and $z$ directions, respectively. Assuming the sample is rigid within a small local region, the OCT image of the displaced sample can be expressed as a convolution of a nondisplaced sample structure and displaced PSF as $s(\mathbf{r}+\mathrm{\Delta }\mathbf{r})=\eta (\mathbf{r})*h(\mathbf{r}+\mathrm{\Delta }\mathbf{r})$, where $\eta (\mathbf{r})$ is the sample structure, $*$ represents convolution, and $\mathrm{\Delta }\mathbf{r}(\mathbf{r})$ is written as $\mathrm{\Delta }\mathbf{r}$ for simplicity. Therefore, the amplitude of the complex correlation coefficient (ACCC) between two OCT signals before and after displacement can be expressed as the amplitude of the correlation between $h(\mathbf{r})$ and $h(\mathbf{r}+\mathrm{\Delta }\mathbf{r})$ as $|\rho (\mathbf{r};\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z,w_l,w_z)|=\exp [-(\mathrm{\Delta }{x}^2+\mathrm{\Delta }{y}^2)/w_l^2-\mathrm{\Delta }{z}^2/w_z^2]$. In this Letter, the OCT images before and after the displacement are referred to as the reference and target images, respectively.

The goal of the present method is to determine the local displacement at the OCT plane, ${\mathrm{\Delta }\mathbf{r}(x,y,z)|}_{y=0}$, from the measured ACCCs. Because our mathematical model for an ACCC has five unknown quantities ($\mathrm{\Delta }x$, $\mathrm{\Delta }y$, $\mathrm{\Delta }z$, $w_l$, and $w_z$), five equations are required to determine $\mathrm{\Delta }\mathbf{r}$. The five equations are obtained by computing local ACCCs between a target image and five digitally shifted reference images such that ${\mu }_0(x,z)\equiv |\rho (\mathrm{\Delta }x,\mathrm{\Delta }z)|$, ${\mu }_1(x,z)\equiv |\rho (\mathrm{\Delta }x+\delta x,\mathrm{\Delta }z)|$, ${\mu }_2(x,z)\equiv |\rho (\mathrm{\Delta }x-\delta x,\mathrm{\Delta }z)|$, ${\mu }_3(x,z)\equiv |\rho (\mathrm{\Delta }x,\mathrm{\Delta }z+\delta z)|$, and ${\mu }_4(x,z)\equiv |\rho (\mathrm{\Delta }x,\mathrm{\Delta }z-\delta z)|$, where the first and second parameters of $\rho$ are substituted into $\mathrm{\Delta }x$ and $\mathrm{\Delta }z$, respectively. $\mathrm{\Delta }x(x,z)$ and $\mathrm{\Delta }z(x,z)$ are respectively written as $\mathrm{\Delta }x$ and $\mathrm{\Delta }z$ for simplicity. These equations represent the ACCCs between the target image and reference images with digital image shifts of (0, 0), $(-\delta x,0)$, $(\delta x,0)$, $(0,-\delta z)$, and $(0,\delta z)$. The values of these equations, i.e., the ACCCs at the OCT plane ($y=0$), are obtained from OCT signals, as described later.

By solving this equation system, $\mathrm{\Delta }x(x,z)$, $\mathrm{\Delta }z(x,z)$, $|w_l(x,z)|$, $|w_z(x,z)|$, and $|\mathrm{\Delta }y(x,z)|$ are obtained as

Note that all digitally shifted reference images are computed from the same reference image; hence, only one measurement is required for the reference image. The amount of digital shift should be smaller than the size of the PSF. In our particular implementation, this was one pixel; hence, $\delta x\sim 1.46\mathrm{\mu m}$ and $\delta z=2.6\mathrm{\mu m}$.

The measurement accuracy of the present method is highly dependent on the estimation accuracy of the ACCCs. Hence, noise-effect-corrected ACCC, defined as follows, was used as an implementation of $|\rho |$ [3,7]:

Figure 1 shows the schematic of the OCT system used for the experiments, which is the same as the one described in Ref. [3]. It is a spectral domain OCT (SD-OCT) that uses an SLD light source (S-1020-B-7, Superlum Ltd., Ireland) with a center wavelength of 1.02 μm and a spectral bandwidth of 100 nm. The interferometer is a fiber Michelson interferometer with a $50/50$ fiber coupler. The interference signal is detected by a high-speed spectrometer comprising an InGaAs line sensor (SUI1024LDH2, Sensors Unlimited, Inc., Goodrich, NC) at a line rate of 91,911 A-scans/s. The probe beam power is 1.86 mW on a sample. The sensitivity was measured to be 95.2 dB, and the measured axial resolution was 5.9 μm in tissue.

 

Fig. 1. OCT system used in the experiments: (a) sample arm for phantom measurements and (b) sample arm for porcine eye measurement.

Download Full Size | PDF

Figure 1(a) shows the sample arm configuration used for tissue phantom measurements. The tissue phantom, which was made from gelatin powder (10 g), fat emulsion (20% Intralipid, 1 g), and 50 ml of hot water, was mounted on a three-axis motorized translation stage (Max301/M, Thorlabs, Inc., NJ). The slow axial drift of the OCT signal was corrected by a glass plate that was placed in front of the sample [3]. The focal length of the objective was 16 mm, and the lateral $1/e^2$ spot radius was 10 μm. The objective was selected so that the dioptric power is similar to that of a porcine eye. The imaging area was 1.5 (lateral) $\times 1.3$ (depth) mm, which is $1024\times 512$ pixels.

Figure 1(b) shows the sample arm configuration used for the measurement of a porcine retina under laser coagulation. The porcine eyes were enucleated within 12 hours after death and investigated within 12 hours after the enucleation. A custom-made zero-diopter contact lens was attached to the cornea to prevent dehydration. The coagulation laser (532 nm wavelength, GYC-1000, Nidek Co., Ltd., Aichi, Japan) was installed in the sample arm. The coagulation beam and OCT probe beam were combined using a dichroic mirror. Furthermore, the coagulation beam and the center axis of the OCT probe beam were coaxially aligned. The power of the coagulation laser was 200 mW, and the theoretical beam diameter on the retina was 83 μm.

The maximum measurable displacement was expected to be around half the width of the PSF [3]. Hence, it was expected to be around 10 μm in the lateral direction and about 5 μm in the axial direction.

We evaluated the accuracy and repeatability of the proposed method by repeating the translation-and-measurement six times. The tissue phantom was translated into either in-plane or out-of-plane directions until the total displacement reached 3 μm. The translation was done with a step size of 0.2 μm, and 16 B-scans were obtained at each step. The spatially resolved displacements were computed from an OCT image taken before translation (reference image) and images at each step (target images). The measured displacements were averaged within a 100 (lateral) $\times 50$ (axial) pixel region of interest, as indicated by the red box in Fig. 2(a). The means and standard deviations of the six mean displacements were calculated at each step.

 

Fig. 2. Representative B-scan of (a) the tissue phantom. The yellow arrow indicates a signal from the glass plate. (b) Means and (c) standard deviations of the mean in-plane displacements. The red crosses and the blue boxes indicate the in-plane lateral and axial translations, respectively. (d) Means and (e) standard deviations of the mean out-of-plane displacements.

Download Full Size | PDF

The measured means and standard deviations of the six mean displacements are plotted as a function of the translation stage position in Figs. 2(b)–2(e). The absolute difference of the means of the six mean displacements from the theoretical displacement was less than 0.15 μm for the in-plane lateral, 0.22 μm for the in-plane axial, and 0.32 μm for the out-of-plane directions. The root mean squared errors were fewer than 0.33 μm for in-plane lateral, 0.23 μm for the in-plane axial, and 0.33 μm for the out-of-plane displacements. The standard deviations of the mean displacements among six repetitive measurements were fewer than 0.36 μm for the in-plane lateral, 0.10 μm for the in-plane axial, and 0.19 μm for the out-of-plane directions. As was expected, the measured in-plane displacements were linear to the set displacement. The measured out-of-plane displacement was also linear, except in the region close to zero displacement.

We measured the thermal dynamics of the laser irradiated porcine retina. Multiple retinal OCT images were obtained by sequentially scanning a same location. The displacement was computed between successive B-scan images. Figure 3 shows the dynamic change of the displacement. The coagulation laser irradiated the sample for 20 ms, where the starting time of the irradiation is denoted as $t=0\mathrm{ms}$. The first to third rows represent the in-plane lateral, in-plane axial, and out-of-plane displacement maps, respectively. Note that the in-plane displacements were set to zero if $w_l$ or $w_z$ took imaginary values. Similarly, the out-of-plane displacement was set to zero if $w_l$, $w_z$, or $\mathrm{\Delta }y$ took imaginary values. The fourth row shows the displacement vectors superimposed on the intensity image. The transparency of the vector arrows represents the magnitude of the displacement, and the color of the arrow shows the orientation of the displacement as represented by a color wheel.

 

Fig. 3. In-plane lateral (first row), in-plane axial (second row), and out-of-plane lateral (third row) displacement maps, as well as the displacement vectors superimposed on the intensity image (fourth row) at different time points during coagulation laser irradiation.

Download Full Size | PDF

It was found that the neural retina expanded radially during laser irradiation and was slowly reinstated after irradiation. In contrast, an area of tissue close to the retinal pigment epithelium (RPE) presented large out-of-plane displacement. However, this could be an artifact caused by alterations of the tissue micro-structure. This is suggested by the fact that the out-of-plane displacements did not decrease, even long after the laser irradiation. These findings suggest the following dynamic process: At first, the coagulation laser is mainly absorbed by the RPE, then the tissue micro-structure is altered and expanded. Finally, the expanded tissue pushes the neural retina in the upward direction.

The present method has several advantages over conventional methods, based on correlation coefficients. First, the present method is immune to noise because of the noise-effect-corrected sample correlation coefficients. Second, the method provides not only the in-plane bidirectional displacement, but also the magnitude of the out-of-plane displacement that was ignored in our previous method [3]. Third, the method is faster and computationally more efficient than speckle tracking and phase-sensitive measurement. Speckle tracking requires a cross-correlation calculation for each local region in the image. In addition, sub-pixel displacement measurements by speckle tracking require computationally intense image interpolation to solve an optimization problem. In the proposed method, however, the computation of local correlation coefficients is required only five times. In addition, computationally intense sub-pixel image shifting is not required. The conventional phase-sensitive measurement requires phase unwrapping if the phase shift between two images exceeds $\pm \pi$-rad range, and this frequently occurs in practical cases [2,3]. Conversely, the proposed method does not use phase information, so it is free from phase wrapping. In our particular implementation, the processing time for a single set of five maps, including $\mathrm{\Delta }x(x,z),\mathrm{\Delta }z(x,z),w_l(x,z),w_z(x,z)$ and $|\mathrm{\Delta }y(x,z)|$, is 1.3 s with a LabVIEW program running on a PC with Intel Core i7 980 CPU (six physical cores, 12 logical cores, 12 M L3 Cache, 3.33 GHz clock speed) and 24 GB memory. The main processing routine uses only a single thread. Although this computational time does not allow real time processing, it would be easily accelerated with multi-thread or GPU processing. Therefore, the presented method is suitable for daily clinical use.

The proposed method has some limitations. First, the basement assumptions are violated in some circumstances. For example, alteration of the micro-structure and deformation within the kernel violate the assumption that the sample is rigid within the kernel. Another assumption violation is caused by aberrations, which cause the PSF to be non-Gaussian and violate the assumption of a Gaussian-shaped PSF. These violations could partially account for occasional imaginary values of $w_l$ and $w_z$. Second, in future in vivo applications, involuntary eye motion causes large displacement that exceeds the maximum measurable range of the method. However, such large displacement can be corrected by using a standard image registration algorithm before applying the proposed method. In addition, a laser induced tissue displacement is expected to be fewer than a few micrometers between successive B-scans [2]. Therefore, the in vivo thermal dynamics could be measurable by combining the image registration algorithm with this proposed method.

We developed a method to measure the bidirectional in-plane displacements and the magnitude of out-of-plane displacement using the correlation coefficients of OCT signals. The method measured displacements with an accuracy better than 0.32 μm and repeatability better than 0.36 μm. In addition, it successfully visualized the dynamics of a porcine retina during coagulation laser irradiation.