Two-dimensional micro-displacement measurement for laser coagulation using optical coherence tomography

: To improve the reproducibility of photocoagulation, the ability to quantitatively monitor the thermal change of laser-irradiated retinal tissue is required. Recently, optical coherence tomography has enabled non-invasive and non-contact monitoring of the tissue structural changes during laser irradiation. To further improve the capability of this technique, a method is proposed to measure tissue displacement by simultaneously using Doppler phase shifts and correlation coefﬁcients. The theoretical approach for this method is described, and its performance is experimentally conﬁrmed and evaluated. Finally, lateral and axial displacements in the laser-irradiated retinal tissues of an enucleated porcine eye are observed. The proposed method is found to be useful for further understanding the direct thermal response of laser-irradiated retinal tissue.


Introduction
Photocoagulation is an effective laser treatment that has been successfully used for eye diseases such as diabetic retinopathy [1,2], retinal vein occlusion, [3,4], and retinal detachment. However, there are iatrogenic complications owing to this treatment, with the most significant drawback being excessive laser damage due to overtreatment. To prevent such damage, it is required that laser irradiation parameters be set at a level that is less than the damage threshold with high reproducibility.
In practice, this requirement is challenging because the laser treatment has poor reproducibility. In most cases, optimization of the laser irradiation treatment relies upon the subjective evaluation of the physician, and relying solely upon the basis of the visibility of whitening lesions. In addition, the scattering and absorption properties of the ocular media, including the pigments, varies spatially, temporally and by the individual, and it makes the process more difficult. Because there are many uncertainties concerning tissue alteration during laser treatment, a quantitative monitoring technique is required.
Recently, real-time and non-invasive monitoring methods have shown potential for improving the reproducibility of results. For example, optoacoustic signals are being used to measure the laser-induced temperature increase [5], which has enabled temperature controlled photocoagulation [6,7]. Another example is the use of optical coherence tomography (OCT) [8]. OCT is used to measure not only the depth-resolved scattering changes but also the tissue displacements that occur because of the thermal expansion induced by the laser coagulation with the help of phase-sensitive measurements [9,10].
Some of the advantages that OCT has over optoacoustic measurements are, first, that OCT enables a non-contact measurement, which is important in routine clinical practice. Second, OCT visualizes a spatially-resolved structural change, which is necessary to observe localized tissue alteration. Third, OCT can directly monitor the thermal changes of the tissue including mechanical displacement, where such displacements can be measurable with the help of phase-sensitive measurements [9][10][11][12], speckle tracking [13], digital image correlation [14], and decorrelation [15,16]. Therefore, OCT is a promising tool for thermal expansion measurements.
In this paper, we propose a method to measure two-dimensional (2D) micro-displacements using OCT and demonstrate thermal expansion measurements during laser irradiation. It is expected that this method can measure the nano-scale axial and sub-micron-scale lateral displacements simultaneously and be more sensitive than the conventional speckle tracking [13]. The rest of the paper is organized as follows. First, the method is described theoretically. In this theory, Doppler phase shifts and correlation coefficients are simultaneously used to measure in-plain displacements. Because we can assume that thermal expansion/contraction is axially symmetric around the axis of laser irradiation, these can be monitored by measuring in-plane displacements in the plane containing the axis of laser irradiation. This creates a simple model to describe the lateral and axial displacements. In addition, we propose a calibration method to estimate the resolution parameters using the OCT cross-sectional (B-scan) image, where the resolution parameters are dependent upon the spatial location of the sample. The spatiallyresolved displacements are then measured using this calibration data and the local analysis of the correlation coefficient and the Doppler phase shift. Second, the performance of correlation coefficient-based displacement measurement is experimentally evaluated using a tissue phantom and a translation stage. Finally, the thermal expansion measurement is demonstrated during laser irradiation using the retinal tissue of an enucleated porcine eye.

Theory
By assuming a laser beam with a Gaussian profile on the sample and a Gaussian temporal coherence function of the light source, the point spread function (PSF) of the OCT system can be expressed as where r = (x, y, z) is the position vector in Cartesian coordinates, x and y indicate the lateral positions and z is the depth position, w l is the lateral 1/e 2 spot radius of a Gaussian probe beam, w z is the 1/e 2 width of the axial PSF, z 0 is the zero-delay depth position of the interferometer, and k c is a wave number corresponding to the center wavelength λ c of the light source. The OCT signal of a sample with a rigid-body displacement can be expressed as a convolution of a non-displaced sample structure and a displaced PSF, h(r + Δr), where the displacement vector is given by Δr = (Δx, Δy, Δz). Here the complex correlation coefficient (CCC) between two OCT signals with and without the rigid-body displacement can be expressed as the CCC of the original PSF and the displaced PSF: This equation indicates that the amplitude of the CCC decreases as a function of sample displacement, and it follows Gaussian. In the application of OCT for laser coagulation, we can assume that the thermal expansion/contraction is axially symmetric around the axis of the laser irradiation, as shown in Fig.  1. We also assume that an OCT B-scan is taken in a plane containing the axis of laser irradiation. Using these assumptions, the out-of-plane displacement Δy is negligible in this B-scan plane, and the thermal expansion can be monitored by measuring only the in-plane displacement (Δx, Δz). Thus, the lateral displacement is given by manipulating Eq. (2) such that where |ρ h | is the amplitude of the CCC (ACCC) between OCT images taken before and after the rigid-body displacement. Equation (3) suggests that the lateral displacement can be determined by measuring the ACCC and the axial displacement Δz if the resolution parameters (w l , w z ) are known.
The axial displacement Δz can be obtained by taking the phase of the CCC, Δφ ≡ ρ h , which is equivalent to the Doppler phase shift. The axial displacement is given by an unwrapped Doppler phase shift: where n is the refractive index of the sample, and U [ ] represents a phase unwrapping operator.
In this application, some quantities, including Δx, Δz, ρ h , and Δφ , should be considered as functions of x and z because the amount of displacement varies across space. In addition, the resolution parameters w l and w z should also be treated as functions of x and z because these parameters can be altered by depth-dependent defocus, system and sample aberrations, scanning nonlinearity, and dispersion [15]. Hence, Eqs. (3) and (4) should be modified as Our purpose is to determine the spatially-resolved maps of the displacement that occurs between before and after laser irradiation based on these equations. Here we have assumed the displacement is no longer macroscopically rigid, but is rigid in a small region in which the CCC is determined. With this assumption, the ACCC |ρ h (x, z)| and the Doppler phase shift Δφ (x, z) are obtained from OCT signals within a small region centered at (x, z) that have the size K x (lateral) × K z (axial) pixels, as described in Sections 3.1 and 3.2. In addition, we determine the spatially-resolved resolution parameters w l (x, z) and w z (x, z), as described in Section 3.3. These equations thereby provide 2D displacement maps of a sample.

Correlation coefficient
In this section, we describe how the ACCC |ρ h (x, z)| is obtained from the measured OCT images. It should be noted that the ACCC obtained from measured OCT images is affected by noise, while the |ρ h (x, z)| in Section 2 has been assumed to be noise-free. Therefore, we need to properly model the noise and correct for it. Initially, we theoretically derive a relationship between population correlation coefficients with and without noise, after which we describe a method to estimate the ACCC at each pixel in a pair of OCT images.

Population correlation coefficient under noise
In OCT images, a noise-free signal at a particular pixel s(x, z) can be expressed as the convolution of the PSF h(x, z) and the summation of the complex field reflectivities of randomly distributed scatterers within the pixel η(x, z): where * represents the convolution. If a sufficiently large number of scatterers contribute to the summation, the OCT signal becomes a fully-developed speckle and is considered a stochastic zero-mean complex circular Gaussian variable S [17,18]. Note that the speckle is considered not as noise but as the varying nature of OCT signal. In addition, when additive noise N in the measurement is also a stochastic zero-mean complex circular Gaussian variable and is independent of the OCT signal S, the resulting measured signal, which is given by the sum G = S + N, is also a zero-mean complex circular Gaussian variable. The realization of a measured signal in the noise is denoted as g(x, z) in this paper.
At first, the complex population correlation coefficient between two noise-free OCT signals S 1 and S 2 is defined as where E [ ] represents the expectancy and the superscript * denotes the complex conjugate. The amplitude of ρ S represents the population correlation coefficient of the noise-free OCT signals.
Similarly, the complex population correlation coefficient between the measured signals G 1 = S 1 + N and G 2 = S 2 + N is defined as where the amplitude of ρ G represents the population correlation coefficient of the OCT signals under noise. By defining the expected signal-to-noise ratio (SNR) of each measurement as SNR i = E |S i | 2 /E |N| 2 where i = 1, 2, the relationship between the population correlation coefficients with and without noise is written as With this formulation, we can estimate the population correlation coefficient of the noise-free OCT signals from that of OCT signals under noise. A more detailed derivation and formulation are available in Ref. 19. Note that ρ i (i is an arbitrary subscript) in this paper represents the CCC, while the same symbol represents the ACCC in Ref. 19. The population correlation coefficients described in this section have provided a framework to understand the relationship between the noise-free and under-noise population correlation coefficients. In the following section (Section 3.1.2), we describe a method to experimentally determine the noise-free correlation coefficients.

Sample correlation coefficient under noise
With this method, the spatially-resolved displacement map is determined using the local correlation coefficients computed from measured OCT signals. This complex correlation coefficient to be measured is denoted as a complex sample correlation coefficient ρ {s,g} , while the complex population correlation coefficient is denoted as ρ {S,G} .
To determine the 2D displacement maps, multiple OCT B-scans are taken by scanning the same position on the sample repeatedly. The measured OCT signal of the nth B-scan, which suffers from noise, is described as g n (x, z). The spatially-resolved in-plane displacements are calculated from a reference B-scan g R (x, z) that is selected from the B-scans taken before the displacement, and a target B-scan g T (x, z) that is selected from the B-scans taken after the displacement.
To obtain a spatially-resolved noise-free sample correlation coefficient |ρ s (x, z)| between the reference B-scan g R (x, z) and the target B-scan g T (x, z), we first compute the complex sample correlation coefficient ρ g (x, z) under noise, which is given by where ξ and ζ respectively represent the lateral and axial positions in a kernel in which the local sample correlation coefficient is computed, and ∑ ξ ∑ ζ represents the summation over the kernel.
Because it is assumed that the noise-free and under-noise sample correlation coefficients follow the same relationship between population correlation coefficients (Eq. (10)), the noisefree sample correlation coefficient |ρ s (x, z)| is determined by where SNR R (x, z) and SNR T (x, z) are the SNRs at each pixel in the target and the reference B-scans. The SNRs are given by where n = R or T , and N is the amplitude of the noise floor recorded before the measurement without the backscattered light from the sample. Although the noise floor is dominated only by the optical noise of the reference beam and the detection noise, it rationally expresses the expectancy of the noise amplitude during a sample measurement. It is because the backscattered light from the biological sample is significantly smaller than that of the reference beam, and it does not become a noise source during the measurement.
It should be noted that this formulation is derived by assuming a fully-developed speckle formation, though the measured OCT signals does not always form a fully-developed speckle. Therefore, it is not guaranteed that the relationship between the noise-free and under-noise sample correlation coefficients is preserved, especially if the sample has low scattering from scatterers with similar refractive index [20]. In addition, the SNR n (x, z) determined by Eq. (13) does not represent the true SNR at (x, z). This is because the noise is stochastic variable in practice, while N in Eq. (13) is an expectancy of noise amplitude determined by a root-meansquare of the realizations of the stochastic noise and, therefore, ρ s (x, z) sometimes becomes larger than 1. The departure of SNR n (x, z) from the true SNR becomes large if the pixel has low signal strength. Despite of these limitations, Eq. (12) statistically provides a more rational estimation of the noise-free correlation coefficient than traditional definition of the sample correlation coefficient [Eq. (11)].

Noise-free sample correlation coefficient and ρ h (x, z)
In principle, the noise-free sample correlation coefficient |ρ s (x, z)| can be substituted into the ACCC (|ρ h (x, z)|) in Eq. (5) to determine the lateral displacement. However, in practice, |ρ s (x, z)| is slightly smaller than |ρ h (x, z)|, mainly because the instability of the lateral scanning of the OCT scanner reduces the correlation coefficient and the space-variant scattering property of the sample affects the correlation coefficient. Therefore, the departure of |ρ To correct this departure, the relationship between |ρ s (x, z)| and |ρ h (x, z)| is modeled as is the spatially-resolved correction factor and represents the effects from scanning instability and imperfection of the signal model. This correction factor is determined by using two successive B-scans, g R (x, z) and g R+1 (x, z), acquired before the displacement. Because there is no displacement between these B-scans, |ρ h (x, z)| becomes unity and the correction factor is defined as where ρ s (x, z) R to R+1 is the noise-free sample correlation coefficient between the two successive B-scans g R (x, z) and g R+1 (x, z), and is computed using Eq. (12).
It should be noted that any lateral displacement is discarded from further analysis if the corrected correlation coefficient γ(x, z) |ρ s (x, z)| exceeds 1, which can occur because of an imperfect estimation of SNR, as discussed in Section 3.1.2. In our implementation, γ(x, z) is computed together with the resolution parameters w l (x, z) and w x (x, z) using regression analysis of the correlation coefficients. The details of this process are described in Section 3.3.
Note that this correction factor corrects static degradation of the correlation coefficients, which is assumed to be constant through the series of B-scans taken during the measurement. On the other hand, Eq. (12) corrects for the dynamic departure of ρ s (x, z) from ρ g (x, z) caused by an alteration of the SNR among the B-scans.

Doppler phase shift
The Doppler phase shifts for each pixel Δφ (x, z) can be estimated by taking the phase of the complex sample correlation coefficients ρ g (x, z). Namely, the Doppler phase shift between the reference B-scan g R (x, z) and the target B-scan g T (x, z) is given by Note that this equation is equivalent to the Kasai autocorrelator [21]. This Doppler phase shift is unwrapped and substituted into Eq. (6) to obtain the axial displacement. In our particular implementation, the 2D phase unwrapping is performed using a quality-guided path-following method [22]. The unwrapped phase at the edge of each OCT B-scan image is used as a reference and forced to be zero. After this unwrapping, the positive and negative axial displacements respectively represent upward and downward displacements in the OCT B-scan image.

Resolution parameters
As discussed in Section 2, the resolution parameters w l (x, z) and w z (x, z) are not constants over the image field. This is not only because of optical aberrations in the system and the nonlinearity in the lateral scanning, but also because of optical aberrations induced by the sample. In addition, the scattering property of the sample can affect the resolution parameters. Namely, the simplest correlation coefficient model (Eq. (2)) implicitly assumes a fully-developed speckle, while the measured OCT signal does not always form the fully-developed speckle. This departure from the implicit assumption is particularly significant for a sample with low scattering from scatterers with similar refractive index [20], such as a retina. Therefore, w l (x, z) and w z (x, z) should be determined for each sample and for each set of measurements.
In our method, these resolution parameters are determined by using successive B-scans of g R (x, z) and g R+1 (x, z) which are both taken before the displacement. At first, g R+1 (x, z) is digitally shifted in the lateral and axial directions with sub-pixel resolution using a zero-padding interpolation. At each shift, a correlation map |ρ s (x, z; Δx, Δz) R to R+1 | is obtained between g R (x, z) Calculate correlation coefficients between the reference and digitally shifted B-scan.
Calculate the sample correlation coefficients Obtain the axial displacement The correlation maps are then fitted using a mathematical model: and γ(x, z), w l (x, z) and w z (x, z) are determined for each position in the B-scan at which the fitting is performed. In our implementation, the fitting is performed with a logarithmic version of Eq. (17), where the model equation is treated as a 2D linear equation of Δx 2 and Δz 2 , and the primary fitting parameters are 1 γ(x,z) , w 2 l (x, z), and w 2 z (x, z) rather than γ(x, z), w l (x, z), and w z (x, z).
Finally, w l (x, z) and w z (x, z) are used as resolution parameters and γ(x, z) is used as the correction parameter of Eq. (15).

Processing flow
In this approach, the processing flow consists of two steps; calibration and displacement measurement. The overall processing flow is summarized in a chart shown in Fig. 2.
In the calibration step, we determine the correction parameter γ(x, z) and the resolution parameters w l (x, z) and w z (x, z), as described in Section 3.3. The digital shift used in this calibration is performed with a shifting range of ± 3 μm both for the axial and lateral directions, with a shifting resolution of 0.2 μm. The digital shift is performed only in the lateral and axial directions and not in the diagonal directions. The noise-free sample correlation coefficient at  each pixel is computed using the neighboring pixels in a 7 × 7 pixel kernel. The shifting range is determined based upon the maximum measurable range, to be discussed in Section 3.5. Figure 3 shows an example of the estimated maps of the inverse of the correction parameter 1 γ(x,z) and the resolution parameters w l (x, z) and w z (x, z). The samples used for these maps are a tissue phantom and an ex vivo porcine retina. It is evident that the correction and resolution parameters are not uniform among the fields, which emphasizes the necessity of the calibration. After determining these parameters, the displacement maps between a reference B-scan and each target B-scan are obtained. Here the lateral displacement is obtained by Eq. (15) and the axial displacement is obtained by Eq. (6).
It should be noted that the displacement measurement is performed for each B-scan, while the calibration process is performed only once for each measurement session, and therefore only once for a series of B-scans. Therefore, the computation load of the calibration process is less problematic than that of the displacement measurement.

Maximum measurable range of correlation analysis
To determine the digital shifting range of the calibration process (Section 3.4) and thereby optimize the resolution parameters, we analyze the maximum measurable range of the correlation coefficient-based displacement measurement. For this analysis, a tissue phantom is measured using a 1 μm-spectral domain OCT (SD-OCT) system. The tissue phantom is made from gelatin powder (10 g), fat emulsion (20% Intralipid, 1 g) and 50 ml of hot water. The details of the SD-OCT system are described in Section 3.5.
An example of the OCT B-scans of the tissue phantom is shown in Fig. 4(a). The noisefree sample correlation coefficients |ρ s (x, z) R to R+1 | are computed between this B-scan and its digitally-shifted versions. Four observation points (P1-P4) are selected randomly, as indicated in Fig. 4(a). The noise-free sample correlation coefficients at these points are plotted in Figs. 4(b) (lateral shift) and 4(c) (axial shift).
The results show that the noise-free sample correlation coefficients monotonically decrease as the absolute shift increases if the shift is sufficiently small. However, the noise-free sam-  ple correlation coefficients begin to oscillate if the shift is too large, and it becomes no longer reliable. This oscillation may be an effect of the structural and optical decorrelation, and the resolution parameters should therefore be determined by using the data only within the monotonic region. In our particular system, the monotonic region is typically within ±10 μm for lateral and ±5 μm for axial directions. Therefore, the fitting range is set as ±3 μm for both the lateral and axial directions, which is safely within the monotonic region.

Spectral domain optical coherence tomography
The schematic of the SD-OCT equipped with a coagulation laser is shown in Fig. 5. The SD-OCT system uses a superluminescent diode 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 a InGaAs line sensor (SUI1024LDH2, Sensors Unlimited, Inc., Goodrich, NC) at a line-rate of 91,911 A-scans/s. This interferometer is equivalent to that described in Ref. 23, and the details of its optical design are described in Ref. 24.
The probe beam power is 1.86 mW on a sample, and the sensitivity was measured to be 93 dB. Because the recoupling loss from the sample to the fiber tip of the sample arm was 2.2 dB for the sensitivity measurement, the real sensitivity was estimated to be 95.2 dB. The axial resolution measured by using a mirror reflector was 5.9 μm (full-width at half maximum (FWHM)) in tissue. The equivalent 1/e 2 width of the axial PSF was 5.0 μm (1/e 2 width = FWHM/ √ 2 ln 2). The mean of estimated axial resolution parameters (see Section 3.3) of the tissue phantom was 5.36 um.

Sample arm configuration
Two sample arm configurations were used in our experiments. The first configuration, shown in Fig. 5(a), was used for phantom imaging and the second configuration, shown in Fig. 5(b), was used for ex vivo porcine retinal imaging using irradiation of the coagulation laser.
In Sections 3.5 and 4.1, a tissue phantom was measured using the first configuration [ Fig.  5(a)]. The tissue phantom was mounted on a three-axis motorized translation stage (Max301/M, Thorlabs, Inc., NJ) which has an on-axis accuracy of 1 μm and a maximum translation of 20 μm. The focal length of the objective lens was 16 mm, and the lateral 1/e 2 spot radius of the Gaussian probe beam was 10 μm. The Gaussian beam radius measured by a beam profiler (BP104-IR, Thorlabs, Inc., NJ) was 8.65 μm. The mean of the estimated lateral resolution parameters (see Section 3.3) of the tissue phantom was 8.74 μm. Multiple OCT B-scans were obtained by a laterally-scanned probe beam. The imaging area was 1.5 mm (lateral) × 1.3 mm (depth), which consisted of 1024 pixels × 512 pixels.
In Section 4.2, the second configuration [ Fig. 5(b)] was used to investigate the retina of ex vivo whole porcine eyes. The porcine eyes were enucleated within 12 hours after death and investigated within 12 hours after the enucleation. These eyes were obtained from a local abattoir (Daimon Ltd., Japan). A custom-made zero-diopter contact lens was attached to the cornea to prevent dehydration. Multiple retinal OCT images were measured by sequentially scanning the same location on the retina. Because the dioptric power of the porcine eye is similar to that of the objective lens of the first configuration [ Fig. 5(a)], the same scanning angle of the galvanometric scanner as that of the first configuration provides a similar scanning range. To observe the thermal expansion dynamics of the retinal tissue, the sample arm was equipped with a coagulation laser (532 nm wavelength, GYC-1000, Nidek Co., Ltd., Aichi, Japan). The coagulation laser and the OCT probe beam were combined using a dichroic mirror, which reflects 1 μm-wavelength light and transmits visible light. The coagulation laser and the center axis of the OCT sample beam were coaxially aligned. The power of the coagulation laser is 200 mW, its duration is 20 ms, and the theoretical beam diameter on the retina is 83 μm.

Validation of correlation coefficient-based displacement measurement
In this experiment, we evaluated the performance of the correlation coefficient-based displacement measurement. The tissue phantom was mounted on the translation stage and translated into either the axial or lateral directions until the total displacement reached 3 μm. Specifically, the translation was done with a step size of 0.2 μm, and 16 B-scans were obtained at each step. At each 0.2 μm translation step, the displacements were computed between a reference B-scan and the obtained B-scans, where the tenth B-scan at the null translation state was used as the reference B-scan. The correction and resolution parameters were calculated from the tenth and eleventh B-scan.
Because the purpose of this experiment is to evaluate the correlation coefficient-based displacement measurement, both the lateral and axial displacements were calculated from the noise-free sample correlation coefficients. Namely, for the lateral translation, the displacement was determined by and for the axial translation, the displacement was determined by Equation (18) was obtained by simplifying Eq. (15), and Eq. (19) is obtained by replacing Δx and w l of Eq. (18) with Δz and w z . Note that these equations are valid only when the displacement is restricted to the lateral or axial directions, and are therefore specialized for this validation experiment only. The measured displacements were observed within an arbitrarily-selected region of interest (ROI) sized of 100 pixels (lateral) times 50 pixels (axial). Figure 6(a) shows an example of a B-scan, where the ROI is indicated by a red box. To quantitate the displacement, a mean displacement was calculated for each 0.2 μm-step translation within the ROI, in which any pixels possessing an SNR less than 10 dB and a corrected correlation coefficient γ(x, z) |ρ s (x, z)| larger than 1 were excluded from the analysis.
In addition, this translation-and-measurement series was repeated six times to evaluate the accuracy and repeatability of the results. The means and standard deviations of the six mean displacements were calculated at each 0.2-μm-step translation. Figure 6(b) shows the mean of the mean displacements plotted as a function of the translation stage position. Here, the position of the translation stage was obtained from the position sensor readout provided with the translation stage, and the relationship between the pixel spacing and the readout was calibrated by using the OCT images. The red and blue plots indicate lateral and axial translations, respectively. The gray line indicates the ideal displacement that is equivalent to the position readout. Because the total measurement time of this experiment was long (around 12 s) and the optical path length difference of the interferometer slowly drifted during the measurement, the OCT signal exhibited a significant axial shift of around a few hundred nanometers. This drift was corrected using the calibration signal indicated by a yellow arrow in Fig. 6(a), which was generated by a cover slip placed between the objective and the tissue phantom.
As is expected, the displacements are linear except in the region near zero displacement. This nonlinear region is owing to the rejection of pixels that have a corrected correlation coefficient larger than 1, which causes the measured displacements to be biased and have a positively skewed and peaked distribution at the edges of the measurable range.
In the lateral direction, the measured lateral displacements tend to overestimate translation. The difference of the mean of the mean displacements from the theoretical displacement was 0.22 μm at the 3 μm translation. On the other hand, in the axial direction, the measured axial displacements tend to underestimate the translation. The difference of the mean of the mean displacements from the theoretical displacement was −0.22 μm with the 3 μm translation. These over-and under-estimations could have several causes. First, the accuracy and the stability of the translation stage affects the results. According to the manufacturer's report, the accuracy of the stage is 0.15 μm at the 3 μm translation, and its thermal stability is 1 μm/ • C. Second, the imperfection in our theoretical model also may cause this error. In our model, the PSF was assumed to be Gaussian [Eq. (1)] and the resolution parameters were determined based on this model. However, the PSF can be non-Gaussian in practice, which would cause the error in the displacement measurement. Future modification of the model would improve the measurement accuracy. Third, an imperfection in the coaxiality of the center axis of the OCT sample beam and the translation axis can cause this error. For instance, deviation from the ideal displacement becomes 0.5 μm with the 3 μm displacement if the axis deviation is 10 degrees. Figure 6(c) shows the standard deviation of the mean displacements among six repetitive measurements plotted as a function of the translation stage position. The red and blue plots indicate lateral and axial translations, respectively. For example, the standard deviations of the mean displacement are 0.24 and 0.11 μm for the lateral and axial directions, respectively, with the 3 μm translation. Note that these standard deviations do not only represent the repeatability of the measurement but also the minimum resolvable displacements.
It should be noted that, in our complete measurement algorithm, the axial displacement was determined using the Doppler phase shift. Because it is a phase-sensitive method, it is expected that the accuracy and repeatability are higher than the correlation coefficient-based method.

Thermal expansion of laser-irradiated retinal tissue
The laser-induced thermal changes in an ex vivo porcine retina were investigated. The laser irradiation was applied during the multiple OCT B-scan acquisition. A successive pair of Bscans were selected from the B-scans taken before the laser irradiation. Specifically, the 88th B-scan was used as a reference B-scan, and was used together with the 89th B-scan to determine the correction and resolution parameters. The laser irradiation was started at the 90th B-scan. The time point t is represented relative to the starting time of the laser irradiation (t = 0 ms) and the duration of the laser irradiation was 20 ms.
Figures 7(a)-(f) show a representative B-scan of OCT intensity and its corresponding 2D displacement map, noise-free sample correlation coefficient map, Doppler phase shift map, lateral displacement map, and axial displacement map at t = +24 ms. The 2D displacement map [ Fig. 7(b)] is a pseudo-color image where the color hue represents the directions of the displacement as represented in a color wheel, and the color saturation represents the amount of displacement. The maximum saturation represents a 10-μm displacement and no saturation (gray) represents no displacement. Here, the saturation is forced to be zero if the pixels have a corrected correlation coefficient larger than 1 and an SNR less than 10 dB. The lightness of each pixel represents the OCT intensity. In addition, to enhance the displacement visibility, a maximum filter is applied in which the pixels having the maximum OCT intensity among neighboring 3 × 3 pixels become the new pixel values at the center of the window.
Although the OCT intensity did not show any characteristic change, large Doppler phase shifts were found at the neural retina, as shown in Fig. 7(d). In addition, a large reduction of the correlation coefficients was found at the retinal pigment epithelium (RPE) complex, as shown    in Fig. 7(c). The 2D displacement was clearly visualized by the pseudo-color displacement map, as shown in Fig. 7(b). It is clear that the neural retinal layers at the irradiation site were axially displaced in the upward direction, whereas significant lateral displacement was found at the RPE complex.
The dynamics induced by the laser irradiation are shown in Fig. 8, where the time series of the OCT intensity B-scan, noise-free sample correlation coefficient map, Doppler phase shift map, and 2D displacement map are represented in the first to fourth rows, respectively. As mentioned above, the laser irradiation initiated at t = 0 and lasted for 20 ms. This time series more clearly describes the findings discussed in the previous paragraph. The full sequence of the images is shown in Media 1, in which the duration of the laser irradiation is indicated by the green frame.

Impact of spatially-resolved resolution parameters upon measurement robustness
It is expected that our approach is robust to the effects of optical aberrations, sample structure, scanning non-linearity and dispersion. This is because our method uses spatially-resolved resolution parameters to account for these perturbations.
To confirm this robustness, our approach (method-1) is compared with an approach in which the resolution parameters are predefined constant values among the image field (method-2). The resolution parameters of method-2 are defined based upon the system design parameters, and are taken to be w l = 10 μm and w z = 5.0 μm (FWHM of 5.9 μm).
The single dataset for each lateral and axial displacement obtained in the same experiment described in Section 4.1, are analyzed by these two methods. The displacements are observed at the four randomly-selected positions (P1-P4) indicated in Fig. 9(a). For each 0.2 μm-step translation, nine displacements are computed with nine target B-scans, and the nine displacements are averaged. The standard deviation of the average displacements among these four positions (P1-P4) are also computed as a measure of the robustness.
The average displacements and their standard deviations are plotted in Fig. 9, where the data are plotted as a function of the translation stage position. Figures 9(b) and (c) show the measured lateral displacement using method-1 and method-2, respectively. On the other hand, Figs. 9(d) and (e) show the measured axial displacements using method-1 and method-2, respectively. The red circles indicate the standard deviations, and red crosses (+), green crosses (×), blue boxes, and purple triangles represent the average displacements at P1-P4, respectively. The gray lines represent the ideal displacement that is equivalent to the position readout.
It is evident that method-1 provides a more accurate displacement measure than method-2 both for the lateral and axial directions. In addition, the standard deviation values clearly demonstrate that method-1 has a greater robustness (lower standard deviation) than method-2. Specifically, the standard deviations for both methods are around 0.3 μm in the region close to the null translation. The standard deviations of method-1 remain constant over the whole translation range, while these of method-2 are around two times larger at the 3-μm translation.
These results clearly suggest the greater robustness of our method, and also reinforce the necessity of spatially-resolved resolution parameters.

Time-dependency of spatially-resolved correction factor
The spatially-resolved correction factor γ is modeled as a time-invariant factor, because it is to correct the static degradation of the correlation coefficients due to the scanning instability and imperfection of the signal model.
To validate the time-invariance of γ, the noise-free sample correlation coefficient is computed with different time-intervals between two B-scans without displacement, which is equivalent to  Fig. 10(a) (a tissue phantom) and Fig. 10(c) (a porcine retina), respectively. In the case of the tissue phantom, the 1/γ was around 1.0 and was not drifted for more than 2 s. Although 1/γ of the porcine retina has relatively larger standard deviation, it is still only 0.02 and the 1/γ was also not drifted. Therefore it is rational to treat γ as a time-invariant factor. Figure 11 describes the time course of the lateral and axial displacement of the porcine retina undergoing coagulation laser irradiation. Figures 11(a) and (b) respectively represent the mean displacements at the RPE complex (ROI-1) and at the neural retinal layers (ROI-2), whose ROIs are indicated by the red boxes in Fig. 7(a). The mean displacements are calculated within ROI-1 and ROI-2, where the pixels having an SNR less than 10 dB and having a corrected correlation coefficient larger than 1 are excluded from the analysis. The yellow region in Fig. 11 indicates the duration of the laser irradiation.

Thermal change during laser irradiation
As shown in Fig. 11(b), both the lateral and axial displacements of the neural retina increase during the laser irradiation and decrease at the cooling stage. Although these displacements decrease, it is not perfectly reversible. This reversibility is consistent with the previous report by Müller et al. [10].
However, the lateral displacement at the RPE complex rapidly increases during the laser irradiation and does not decrease after the laser irradiation. These lateral displacements are also significantly larger than the axial displacements. This large and persisting lateral displacement can be explained by two mechanisms. First, the displacement at the RPE complex will be larger     than that of the neural retina because the coagulation laser is primarily absorbed at the RPE. Second, the high laser absorption at the RPE will result in thermal tissue alteration and a consequent tissue deformation (not displacement). Both the displacement and deformation results in a reduction of the correlation coefficient. It should be noted that our theory accounts only for rigid-body in-plane displacements. Therefore, if a deformation and out-of-plane displacement exists, the further reduction of the correlation coefficient causes an overestimation of the lateral displacement. The authors are currently working on a new version of this method, which can distinguish the in-plane displacement from other effects, such as deformations and out-of-plane displacements.

Limitations
The proposed method still has some limitations. First of all, the lateral displacement is overestimated if the sample suffers from deformations and out-of-plane displacements, as partially discussed at the end of Section 5.3. This is an issue for future development.
A second limitation is the ambiguity in the right-left directions of the lateral displacement. Because Eq. (15) provides only the absolute value of the lateral displacement, the direction of the lateral displacement cannot be defined. To obtain the direction, this method can be used in conjunction with a speckle tracking algorithm [13] or solving a minimization problem [25]. Since these extensions of the method would result in a penalty of computation time, their adoption would depend upon the application.
Third, the correlation coefficient reduces with time due to general fluctuations, such as sample deformation [14]. This limits the measurable time separation between the reference and target B-scans. In our laser irradiation experiment, we have limited the time interval between the B-scans to be less than 250 ms, where the reduction of the correlation coefficient is experimentally confirmed to be negligible. Because the pulse duration of the coagulation laser is 20 ms, this measurable time would be sufficient.
Fourth, motion artifacts are crucial in in vivo measurements. Involuntary eye motions cause a reduction of the correlation coefficients. Similarly, a large axial motion and displacement causes phase wrapping in the Doppler measurement. Although a phase unwrapping algorithm can determine the phase offsets (wrapping numbers among pixels), it is still yielding relative phase offsets. To determine the absolute Doppler phase, we need a reference point at which it is sure that no axial motion and displacement exist. In our ex vivo experiment, we use the tissue at the edge of the image as a reference point. To apply this method to in vivo cases, however, a high-performance motion correction algorithm will be required.
Finally, the measurable range is limited, as described in Section 3.5, though this is not a significant problem when the measured displacement is only a few micrometers. For instance, the laser-induced tissue displacement is expected to be less than a few micrometers between B-scans [10]. A future modification containing a protocol that would measure the displacement between successive B-scans, would overcome this limitation.

Conclusions
We have demonstrated a spatially-resolved, rigid-body, in-plane displacement measurement method using the correlation coefficient and the Doppler phase shift. By estimating the correction and resolution parameters in the calibration process, the lateral displacement measurement becomes robust to the optical aberrations, sample structure, scanning non-linearity and dispersion. In our system, the correlation coefficient-based method achieved an accuracy of 0.22 μm and repeatability of 0.24 μm for the lateral displacement of 3 μm. The measurable ranges were less than 10 and 5 μm for the lateral and axial directions, respectively. The proposed method successfully visualizes tissue displacements of laser-irradiated porcine retina and provides use-