1. 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–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 spatially-resolved 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.

2. 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

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

 

Fig. 1 (a) The biological sample is modeled as a cylindrical shape. The coagulation laser is applied at a particular position on the sample. (b) The top cross-sectional view of the sample. The thermal expansion/contraction is modeled as axially-symmetric radial displacements in the x–y plane. (c) The side cross-sectional view of the sample. The thermal expansion/contraction is modeled as in-plane displacements in the x–z plane. The black arrows indicate the displacement vectors. The pink ellipse indicates the main absorber of the sample.

Download Full Size | PDF

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:

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.

3. Implementations

3.1. 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.

3.1.1. 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):

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₁ and S₂ is defined as

Similarly, the complex population correlation coefficient between the measured signals G₁ = S₁ + N and G₂ = S₂ + N is defined as

By defining the expected signal-to-noise ratio (SNR) of each measurement as SNR_i = E [|S_i|²]/E[|N|²] where i = 1, 2, the relationship between the population correlation coefficients with and without noise is written as

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.

3.1.2. 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

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 noise-free sample correlation coefficient |ρ′_s(x, z)| is determined by

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-mean-square 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)].

3.1.3. 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 |ρ′_s(x, z)| from |ρ_h(x, z)| should be corrected.

To correct this departure, the relationship between |ρ′_s(x, z)| and |ρ_h(x, z)| is modeled as |ρ_h(x, z)| = γ(x, z)|ρ′_s(x, z)|, where γ(x, z) 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

Finally, γ(x, z)|ρ′_s(x, z)| is substituted into |ρ_h(x, z)| of Eq. (5) to determine the lateral displacement, such that

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.

3.2. 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

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.

3.3. 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) and g_R+1(x + Δx, z+ Δz), where Δx and Δz respectively represent the amounts of digital shift in the x and z directions.

The correlation maps are then fitted using a mathematical model:

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).

3.4. 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.

 

Fig. 2 The flowchart of the processing steps.

Download Full Size | PDF

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 $\frac{1}{\gamma (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.

 

Fig. 3 The inverse of the estimated correction parameters (a) $\frac{1}{\gamma (x,z)}$ and the estimated resolution parameters (b) ŵ_l(x, z) and (c) ŵ_z(x, z) of the tissue phantom. These of the porcine retina are shown in (d), (e) and (f), respectively.

Download Full Size | PDF

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.

3.5. 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 noise-free 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).

 

Fig. 4 (a) Representative B-scan image of tissue phantom. The red, green, blue and purple crosses indicate the randomly selected positions (P1–P4). (b) The noise-free sample correlation coefficient |ρ′_s(Δx)_{R to R+1}| is plotted as a function of the lateral shift. (b) The noise-free sample correlation coefficient |ρ′_s(Δz)_{R to R+1}| is plotted as a function of the axial shift.

Download Full Size | PDF

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 sample 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.

3.6. 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.

 

Fig. 5 Experimental setup. (a) The sample arm for experimental verification. A cover slip is used to compensate for the drift of the interferometer in Section.4.1. (b) The sample arm for ex vivo experiments.

Download Full Size | PDF

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² width of the axial PSF was 5.0 μm ( $1/e^2\text{width}=\text{FWHM}/\sqrt{2\text{ln}2}$). The mean of estimated axial resolution parameters (see Section 3.3) of the tissue phantom was 5.36 um.

3.7. 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² 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.

4. Results

4.1. 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

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.

 

Fig. 6 (a) Representative B-scan image of the tissue phantom. The red box indicates the selected region of interest. A yellow arrow indicates the calibration signal generated by a cover slip. (b) The means of the mean displacements are plotted as a function of the translation stage position. Red and blue crosses indicate the lateral and axial translations, respectively. (c) The standard deviations of the mean displacements among six repetitive measurements are plotted as a function of the translation stage position. Red and blue crosses indicate lateral and axial translations, respectively.

Download Full Size | PDF

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.

4.2. 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 B-scans 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.

 

Fig. 7 (a) OCT B-scan intensity image at t = 24 ms. Red boxes indicate the regions of interest. Each sub-figure shows (b) a 2D displacement map where the color hue represents the directions of the displacement as represented in a color wheel, (c) a noise-free sample correlation coefficient map, (d) the Doppler phase shift prior to phase unwrapping, (e) the measured lateral displacement, and (f) the measured axial displacement.

Download Full Size | PDF

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.

 

Fig. 8 The OCT intensity (first row), the noise-free sample correlation coefficients (second row), the Doppler phase shifts (third row) and the displacement map (fourth row) are shown at different time points. The full sequence of images is shown in Media 1, in which the duration of the laser irradiation is indicated by the green frame.

Download Full Size | PDF

5. Discussion

5.1. 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.

 

Fig. 9 (a) Representative OCT B-scan image of the tissue phantom. The red cross (+), green cross (×), blue box, and purple triangle indicate the randomly-selected positions (P1–P4). The lateral displacements using (b) method-1 and (c) method-2. The axial displacements using (d) method-1 and (e) method-2. The red crosses (+), green crosses (×), blue boxes, and purple triangles represent the average displacement at P1–P4, respectively. The red circles represent the standard deviations of the measured displacement among P1–P4.

Download Full Size | PDF

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.

5.2. 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 $\frac{1}{\gamma (x,z)}$, and plotted as a function of the time-interval as shown in Figs. 10(b) and 10(d).

 

Fig. 10 The inverse of the estimated correction parameters (a) $\frac{1}{\gamma (x,z)}$ of the tissue phantom and 1/γ at the position indicated by a red cross is plotted as a function of time-interval between two B-scans (b). These of the porcine retina are shown in (c) and (d), respectively.

Download Full Size | PDF

Figures 10(b) and (d) show the time-dependency of 1/γ at a particular position indicated by red crosses in B-scans of 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.

5.3. Thermal change during laser irradiation

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.

 

Fig. 11 The measured lateral and axial displacements are plotted as a function of B-scan ID. Displacements are observed at the (a) RPE complex and (b) neural retinal layer. The yellow region indicates the duration of the laser irradiation.

Download Full Size | PDF

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.

5.4. 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.

6. 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 useful information for the further understanding of the thermal changes that occur during laser coagulation.