Real-time co-localized OCT surveillance of laser therapy using motion corrected speckle decorrelation

Raphaël Maltais-Tariant; Raphaël Maltais-Tariant; Caroline Boudoux; Caroline Boudoux; Néstor Uribe-Patarroyo; Néstor Uribe-Patarroyo

doi:10.1364/BOE.385654

1. Introduction

Laser therapy allows for targeted removal or thermal coagulation of tissues through light absorption by various absorbers in the tissue, such as water, hemoglobin and other macromolecules [1–3]. It is used in the treatment of epithelial cancerous and pre-cancerous lesions such as Barrett’s esophagus [4–6], in eye diseases such as diabetic retinopathy [7–10] and retinal vein occlusion [11,12] and in respiratory papillomatosis, a viral-induced disease in the upper airway [13,14]. However, accurate monitoring is essential to avoid over- or under-treatment. Under-treatment requires repeat visits with the associated cost, while overtreatment generally causes complications requiring surgical treatment [15–17]. Such precise monitoring requires real-time imaging techniques capable of resolving the thin epithelial layers (generally a few hundred microns) under treatment. Furthermore, when the target zone is small and presents mostly subsurface features for guidance, properly targeting the therapy can prove challenging unless treatment is guided by cross-sectional imaging.

Current monitoring methods for laser therapy are limited. Ultrasound imaging has been used to monitor photothermal therapy in cancer treatment, but lacks sufficient resolution to accurately assess treatment depth in the epithelial layers of tissue [18]. Photoacoustic thermography [19] also has limited resolution, requires the injection of nanoparticles, and only provides an indirect measurement of the tissue coagulation. Thermal magnetic resonance imaging has been demonstrated for treatment monitoring, but this approach is costly and presents challenges such as motion artifacts and complex procedures for choosing the optimal monitoring parameters [20–22].

Alternatively, optical coherence tomography (OCT) is a cross-sectional imaging technique with a penetration depth of 1–2 mm and a depth resolution of 1–15 µm [23]. It has been used in identifying pathologies in the epithelium layer of different tissues, including some treatable with laser therapy, such as Barrett’s esophagus [24–26]. Tissue coagulation produces changes in the acquired OCT tomograms which enable the use of OCT to monitor the progression of the coagulation process and determine treatment depth [27–33].

Several methods have been proposed to guide treatment with OCT. Vasculature observed by OCT was used to determine the effect of treatment, but is used to study long-term effects and not real-time monitoring [34]. Changes in tissues scattering properties have been exploited, but it tends to underestimate the therapy depth. The increase in scattering in the superficial layers due protein denaturation causes enhanced attenuation of the light once it reaches deeper layers —casting a “shadow”— which results in an overall decrease in the signal deeper into tissue, inhibiting accurate assessment of the coagulation state of deeper layers [28,29]. Polarization-sensitive OCT can assess the thermal damage directly, but its contrast is limited to birefringent tissues at baseline —tissues that are naturally birefringent—, such as muscles, skin and tendons [35–38]. Doppler OCT combined with analysis of the complex correlation coefficient has also been investigated as a way to measure tissue displacement during coagulation [27,39], however no determination of the level of coagulation of tissue was attempted. An OCT angiography technique known as complex differential variance (CDV) [40] has been employed with success to differentiate coagulated tissues in laser and radio-frequency therapy [29,30]. CDV makes use of the complex-amplitude OCT signal, which adds significant complexity to a real-time implementation of the signal processing. Indeed, many clinical systems lack phase stability: the presence of patient motion or the use of endoscopic probes with scanning inaccuracies such as non-uniform rotation distortion make it necessary to implement several steps of phase jitter correction and motion and distortion correction difficult to achieve in real time [30]. To the best of our knowledge, OCT real-time monitoring capable of dynamically controlling the therapy laser when the target treatment depth has been achieved has not been demonstrated.

In this work, we present a system capable of real-time monitoring of the coagulation depth in epithelial tissue compatible with catheter imaging. A double-clad fiber coupler (DCFC) is used to combine both the OCT and therapy laser in a single double-clad fiber (DCF) for co-registered imaging and therapy [41–44]. Such couplers have been employed with success in combining OCT imaging with other modalities [45–48]. They allow for multimodal single-fiber catheter. This represents an improvement from Ref. [29], where the two lasers were combined using a dichroic mirror, a configuration appropriate for benchtop validation but impractical for clinical use and incompatible with endoscopic systems. We achieve monitoring using the intensity correlation of the speckle, which lowers the signal processing complexity significantly compared to the complex-amplitude speckle processing required with CDV. Furthermore, novel noise- and motion-correction methods are added to simultaneously extend the maximum monitoring depth and improve system robustness to bulk tissue motion. We implemented these signal processing methods in both post-processing and real time. We demonstrate the use of the developed system in a benchtop configuration for real-time coagulation depth monitoring and targeting, where the system stops the therapy laser automatically when the desired treatment depth has been reached, despite the presence of bulk tissue motion. We validated our approach in ex vivo rat tongue and abdominal muscle tissue, with histology confirmation. Our approach opens the possibility for future real-time image-guided laser therapy with high spatial and temporal resolution through a single fiber endoscope.

2. Theoretical model

2.1 Speckle decorrelation during coagulation

Speckle results from the interference of the coherent light reflected or back-scattered off the scatterers present within a turbid sample. In OCT, a given sub-resolution spatial distribution of scatterers with its associated scatterer size and refractive index yields a specific speckle pattern provided there are no intrinsic —such as diffusion— or extrinsic —such as bulk tissue motion— mechanisms affecting their configuration state [49]. However, as tissue undergoes photocoagulation, proteins denaturate, producing a change in the refractive index distribution and, therefore, a change in the coherent sum of individual back-scattered contributions resulting in a temporal evolution of the speckle pattern seen in the tomogram [31,32,39,50]. Lo et al. observed that these variations slow down once the tissue has fully coagulated [29]. By monitoring the speckle pattern dynamics as it transitions from static (before coagulation), to dynamic (during active coagulation) and back to static (after coagulation), it becomes possible to determine the state of each part of the tissue during therapy [29]. The present section will explain our theoretical model. For a direct application, Fig. 2, at the end of this section, shows a flow chart of the formulas used in our approach.

The quantitative analysis of the speckle fluctuations can be performed by determining the temporal correlation of the OCT signal. When relying on the intensity of the signal, the second-order autocorrelation function $g^{(2)}$ is generally used to quantify the correlation. In what follows, we describe how $g^{(2)}$ can be used to create a map of the OCT signal correlation across the field of view. We define a discrete coordinate system $(z,x,y)$ defining the axes in depth, in-plane lateral location, and out-of-plane lateral location, respectively. In the following, we do not consider the out-of-plane dependency, unless otherwise specified.

The second-order autocorrelation function $g^{(2)}$ is defined as the time-delayed correlation of the intensity of an ensemble at a time difference $\tau$. Considering that an OCT system determines the intensity of backscattered light as a function of depth $z$ (first index), lateral position $x$ (second index) and time $t$ (third index), the inter–B-scan autocorrelation function can be written as [51,52]

(1)$$g^{(2)}(z, x, \tau) = \frac{\left\langle \sum_{t = 0}^{\Delta T - \tau} I(k,l,t) \odot I(k,l,t+\tau) \right\rangle _{kl}}{\sum_{t = 0}^{\Delta T - \tau} \left\langle I(k,l,t) \right\rangle \left\langle I(k,l,t+\tau) \right\rangle _{kl}},$$

where $\Delta T$ is the time lapse (measured in number of B-scans) used to calculate the correlation, $\odot$ is the Hadamard product, and $\left\langle \cdots \right\rangle _{kl}$ represents an ensemble average of $I(k,l,t)$ in the $k$ and $l$ indices centered at ($z$, $x$). The OCT system is assumed configured to acquire B-scans as a function of time at the same out-of-plane location; given that Eq. (1) is an inter–B-scan metric, time $t$ is equivalent to the B-scan index. A precise definition of the ensemble will be given in Eq. (3). This is a direct translation of the estimation of the autocorrelation function for a stationary process. However, as explained above, during tissue coagulation the tissue transitions through different states and, therefore, the acquired signal corresponds to a dynamic process. In a dynamic process, $g^{(2)}$ depends not only on $\tau$, but also on $t$. In order to adapt Eq. (1) for a dynamic process to detect the state of coagulation, we limit the temporal extent of the summation to allow us to calculate a time-dependent correlation —assuming the sample is in the same state inside a limited time interval. If we define $\Delta t$ as the number of B-scans used for the instantaneous, estimation of the autocorrelation at time $t$, Eq. (1) becomes

(2)$$g^{(2)}(z, x, t, \tau) = \frac{\left\langle \sum_{m = t}^{\Delta t} I(k,l,m) \odot I(k,l,m-\tau) \right\rangle _{kl}}{\sum_{m = t}^{\Delta t} \left\langle I(k,l,m) \right\rangle \left\langle I(k,l,m-\tau) \right\rangle _{kl}}.$$

Notice that we use $m - \tau$ for our second ensemble since, for the calculation of the correlation in real time, we only have access to the B-scans preceding $t$. To simplify the notation in what follows, we drop the explicit summation and change the definition of ensemble to

(3)$$I_{ens}(z,x,t) \equiv \{I(k,l,m)\},\ \mbox{for} \left\{ \begin{array}{lcl} \forall k & \in & \{z,z +1,\ldots,z + \Delta z\} \\ \forall l & \in & \{x, x+1,\ldots,x + \Delta x\} \\ \forall m & \in & \{t - \Delta t,t - \Delta t+1,\ldots,t\} \end{array} \right.,$$

where $I(k,l,m)$ is the intensity of the pixel at depth $k$, lateral position $l$ and time $m$; and $\{\cdots \}$ denotes the concatenation of the argument into an ensemble. The ensemble $I_{ens}(z,x,t)$ encompasses all pixels at depths $z$ to $z + \Delta z$, lateral positions $x$ to $x + \Delta x$ and times $t$ to $t - \Delta t$. Our time-dependent autocorrelation function becomes

(4)$$g^{(2)}(z,x,t,\tau)=\frac{\left\langle I_{ens}(z,x,t-\tau)\odot I_{ens}(z,x,t) \right\rangle }{\left\langle I_{ens}(z,x,t-\tau) \right\rangle \left\langle I_{ens}(z,x,t) \right\rangle },$$

where the bracket $\left\langle \cdots \right\rangle $ denotes an ensemble average implemented as $\sum _{klm}$. Decorrelated speckle patterns will have $g^{(2)} \rightarrow 1$ while correlated patterns will have $g^{(2)} \rightarrow 2$ for speckle with unity contrast, as is the case of OCT [52]. These theoretical values for $g^{(2)}$ are reached for infinitely large ensembles. Since the ensembles considered here are small, we define a normalized second-order autocorrelation function $g^{(2)}_n$ that accounts for non-unity contrast as

(5)$$g^{(2)}_n(z,x,t,\tau)=\frac{g^{(2)}(z,x,t,\tau)-1}{g^{(2)}_{0}(z,x,t,\tau)-1} + 1.$$

The zero-delay autocorrelation used for normalization $g^{(2)}_{0}$ takes into account the contrast of the speckle in all the samples used in the calculation of the time-dependent autocorrelation function and is thus a function of $\tau$ and $t$ given by

(6)$$g^{(2)}_{0}(z,x,t,\tau) = \frac{\left\langle \{I_{ens}(z,x,t-\tau),I_{ens}(z,x,t)\} \odot \{I_{ens}(z,x,t-\tau),I_{ens}(z,x,t)\} \right\rangle }{\left\langle \{I_{ens}(z,x,t-\tau),I_{ens}(z,x,t)\} \right\rangle ^2}.$$

The term $\{I_{ens}(z,x,t-\tau ),I_{ens}(z,x,t)\}$ is a concatenation of the two ensembles $I_{ens}(z,x,t-\tau )$ and $I_{ens}(z,x,t)$. We ignore any dependency of the calculated value of the autocorrelation function with the ensemble size as we do not compare results among different ensemble sizes [52]. During therapy, we calculate $g^{(2)}_n$ for all positions in the OCT tomogram. Figure 1 shows an example for a given value of $\tau$ and $t$. We elaborate on the choice of $\tau$ in Sec. 4.2. As for $\Delta z$, $\Delta x$ and $\Delta t$, a larger ensemble provides a better estimation of the correlation at the cost of resolution in $z$, $x$ and $t$.

Fig. 1. Overview of the $g^{(2)}$ calculation. Each pixel of the $g^{(2)}$ image (right) is calculated from two ensembles of intensity data points, each ensemble comprising several pixels in $z$, $x$ and among $\Delta t$ B-scans in time. The two ensembles are delayed in time by $\tau$.

Download Full Size | PDF

Fig. 2. Flow chart used to calculate the noise- and motion-corrected second-order autocorrelation function $g^{(2)}_{nc\&mc}(z,x,t,\tau ).$

Download Full Size | PDF

2.2 Noise correction

Noise, being dynamic in time, will lower the correlation coefficient calculated [53]. This effect becomes significant during coagulation when the increased scattering causes stronger attenuation of the signal coming from deeper tissue layers, decreasing SNR. This limits the maximum possible monitoring depth of the therapy [28,29]. As such, it is crucial to implement noise correction methods that account for these effects. We use a novel analytic correction for the noise-induced decorrelation, as proposed by Uribe-Patarroyo et al. [52], which does not suffer from the loss in temporal resolution of previous corrections [53]. Briefly, the shot noise component in the OCT tomogram can be modeled as an additive white-noise contribution to the complex-valued tomogram [54]. The noisy complex-valued signal $G$ is defined as $G = S + N$ where $S$ is the noiseless signal and $N$ is the noise. The noisy signal intensity $I$ is therefore $I = |S + N|^2$. We also define a noiseless signal intensity $T = |S|^2$. Temporarily replacing the indices $(z,x,t-\tau )$ and $(z,x,t)$ with $_1$ and $_2$, respectively, it is possible to write the second-order autocorrelation function of the noisy signal as

(7)$$ g^{(2)}(z,x,t,\tau) = \frac{g^{(2)}_{c}(z,x,t,\tau) + G}{1 + G},\quad g^{(2)}_{0}(z,x,t,\tau) = \frac{g^{(2)}_{c,0}(z,x,t,\tau) + 2G_0}{1 + G_0},$$

where

(8)$$g^{(2)}_{c}(z,x,t,\tau) = \frac{\left\langle T_1 T_2 \right\rangle }{\left\langle T_1 \right\rangle \left\langle T_2 \right\rangle }$$

is the second-order autocorrelation function in the absence of noise, and $g^{(2)}_{c,0}$ is defined in a similar fashion. The parameter $G$ is given by the average SNR of the ensembles, written as $R$:

(9)$$ G(z,x,t,\tau) = \frac{1+R_1+R_2}{R_1 R_2},\quad G_{0}(z,x,t,\tau) = \frac{1+2 R_{12}}{R^2_{12}},$$

where

(10)$$R_j = \frac{\left\langle I_j \right\rangle }{\left\langle |N|^2 \right\rangle } - 1,$$

and $R_{12}$ denotes the SNR of the two ensembles $\{I_1,I_2\}$. It is straightforward to express $g^{(2)}_{c}$ [$g^{(2)}_{c,0}$] as a function of $g^{(2)}$ [$g^{(2)}_{0}$] and obtain noise-corrected versions of Eq. 4. The normalized and SNR-corrected autocorrelation function becomes

(11)$$g^{(2)}_{nc}(z,x,t,\tau)=2\frac{g^{(2)}(z,x,t,\tau)[1+G]-G}{g^{(2)}_{0}(z,x,t,\tau)[1+G_{0}]-2 G_{0}}.$$

The DCF is known to have some level of cross-talk between the single-mode core and inner cladding, which results in ghost images and/or increased noise floor. For this reason, we decided not to rely on the noise floor determined by blocking the sample arm [41]. We thus developed a calibration procedure to accurately determine the noise level, using a tomogram acquired prior to the start of the coagulation process. An additional advantage of this method is the possibility of estimating the SNR for data for which the noise floor is unknown. First, we assumed that the SNRs of both ensembles are identical, $R_1 = R_2$, which proved to be an excellent approximation. This also implied that $G=G_{0}$. We also assumed that the noise intensity remains constant for the remainder of the coagulation, and that the SNR will, therefore, only change due to variations in the signal intensity

(12)$$R(z,x,t) = \frac{\left\langle I_{ens}(z,x,t) \right\rangle }{\left\langle |N|^2 \right\rangle }-1,$$

where $I_{ens}(z,x,t)$ is the $I$ ensemble located at $z$, $x$ and $t$, and it is made explicit that $N$ is independent of $z$, $x$ and $t$. In order to confirm the linear relationship between the $R$ and $\left\langle I_{ens} \right\rangle $, $R$ was estimated for varying $\left\langle I_{ens} \right\rangle $ values. More precisely, we estimated the average SNR of the $j$-th ensemble, $\bar {R}_j$, as a function of its average signal intensity, $\left\langle I_j \right\rangle $, for all the ensembles in the field of view. This calibration was found by solving the optimization problem outlined in Eq. (13). For a static sample with no ongoing treatment, the only source of decorrelation will be the noise. We can therefore expect $g^{(2)}_{nc}(\tau > 0) = 2$ and write

(13)$$\mathop {\arg\,\min}\limits_{_{{{\bar{R}}_j} \in {{\mathbf R}_{ \ge 0}}}} {\left|{g^{(2)}_{nc}(\bar{R}_j, \tau > 0)-2}\right|} \textrm{ for } j = 1,\ldots,m,$$

$m$ being the total number of ensembles. Since the estimation of $\bar {R}_j$ as a function of $\left\langle I_j \right\rangle $ for a single ensemble tends to be inaccurate, ensembles were grouped in $M$ bins: ensembles with similar average intensity were grouped in the $j$-th bin with average intensity $I_j$. Considering the full dynamic range of the signal, we equally distributed all ensembles into $M=100$ bins [thus $j=1\ldots 100$ in Eq. (13)], and and the optimization problem was solved for each bin. We selected the number of bins by finding a good compromise between an accurate estimation of the SNR and a good resolution on the $\bar {R}_j$ vs. $\left\langle I_j \right\rangle $ curve and the time it took to solve all of the optimizations.

The optimization problem itself was solved using a similar approach to the golden-section search method [55]. The function is evaluated at three points, $x_1$, $x_2$, and $x_3$ separated by intervals of $\Delta x$. Step one consists of finding the minimum value between the function evaluated at these three points. In step two, we created a new interval of three points separated by $\Delta x/2$ centered around the minimum calculated in step one. The process is iterated until the interval separating the three points is considered small enough. In our case, this value was 1% of $\left\langle I_j \right\rangle $. This method was chosen because it was found to converge rapidly with our data; however, in principle, any optimization method can be used as the search space seems to have a single global minimum. The curve is used as a lookup table to calculate $g^{(2)}_{nc}(z,x,t,\tau )$ using the proper SNR value, searched for making use of only the current ensemble intensity at each $t$.

2.3 Motion correction

It is well established that tissue bulk motion will affect the value of $g^{(2)}$ as well, thereby leading to inaccurate measurements of the coagulation-induced decorrelation. To address this impact a second correction, based on a lookup table, is applied after the noise correction. The two-dimensional table gives the $g^{(2)}_{nc}$ value caused by motion-induced decorrelation on the horizontal axis and by coagulation-induced decorrelation on the vertical axis. Points lying in between the two axes correspond to decorrelation caused by a combination of both coagulation and motion.

The method works as follows. A look-up table $g^{(2)}_{nc}(k,l)$ is created, in which the row index, $k$, corresponds to the degree of coagulation-induced decorrelation, and the column index, $l$, corresponds to the motion-induced decorrelation. The first row $g^{(2)}_{nc}(k=0,l)$ corresponds to the expected value of $g^{(2)}_{nc}$ for untreated tissue and all possible degrees of tissue motion (given by $l$). For a given experiment, motion is quantified by dynamically calculating $g^{(2)}_{nc}$ in an untreated section of the tissue, which we denote the reference autocorrelation, $\tilde {g}^{(2)}_{nc}$. The value of $\tilde {g}^{(2)}_{nc}$ will be searched for in the first row of the table. The associated column index, $l_{R}$, will define the motion-induced decorrelation for all ensembles in the current frame. This is represented as step 1 in Fig. 3. In step 2, to correct a given $g^{(2)}_{nc}(x)$, the index $k_m$ having the closest match in the column given by $g^{(2)}_{nc}(k,l_R)$ is found, such that $g^{(2)}_{nc}(x)\approx g^{(2)}_{nc}(k_m,l_R)$. In step 3, the noise- and motion-corrected $g^{(2)}_{nc\&mc}(x)$ is found as $g^{(2)}_{nc\&mc}(x)=g^{(2)}_{nc}(k_m,l=0)$, since the values in the first column corresponds to the decorrelation induced by coagulation only.

Fig. 3. Lookup table used for motion correction created from the simulation of coagulation- and motion-induced decorrelation. The vertical axis represents coagulation-induced decorrelation while the horizontal axis represents motion-induced decorrelation. The correction method is as follows: $\tilde {g}^{(2)}_{nc}$ is calculated in an untreated region of the tissue to quantify motion. In step 1. we search for $\tilde {g}^{(2)}_{nc}$ in the first row of the table to determine the motion column index. In step 2. we find the $g^{(2)}_{nc}(x)$ to be corrected for motion in that column and finding its row index. In step 3. we find the $g^{(2)}_{nc\&mc}(x)$ that corresponds to the first value in that row.

Download Full Size | PDF

Since the method uses the correlation to quantify the motion in real time during the therapy, it can also correct for out-of-plane motion. The lookup table in Fig. 3 was created by simulating a noiseless OCT signal exhibiting speckle obtained by randomly moving scatterers —mimicking the coagulation process— along the $l$ index, and uniformly moving scatterers —mimicking bulk motion— along the $k$ index, based on the model of a static sample described by Zaitsev et al. [56]. The method is detailed in Sec. 7. The advantage of simulated speckle is that the motion can be accurately modeled and a considerable quantity of speckle patterns can be simulated, allowing for larger ensembles and thus a more accurate calculation of the correlation.

3. Materials and methods

3.1 Experimental setup

We used a commercial OCT system based on a MEMS-VCSEL wavelength-swept laser at 1310 nm (Thorlabs SL1310V1, Newton, USA) with its imaging module (Thorlabs OCS1310V1, Newton, USA) using a single detector (no polarization-diverse detection) for imaging. We used an Alazar digitizer model ATS9350 to acquire the electronic signal from the imaging module. The laser, with a wavelength range of >100 nm (at -20 dB), has a 16 µm theoretical axial resolution —we determined a 20 µm axial resolution in air experimentally— with an A-line rate of 100 kHz. A CW 532 nm laser (Verdi 10W, Santa Clara, USA) was used as the therapy laser source. A DCFC (Castor Optics DC1300LQ1, Montreal, Canada) was used to couple the therapy laser into the OCT sample arm, which used a single 2.18 m-long DCF. The DCFC enables the delivery of both imaging and therapy through a single fiber, guaranteeing co-registration and greatly simplifying the implementation of our technique in endoscopic probes. The DCF core diameter, cladding diameter, and effective refractive index are all chosen as to match the guided mode waist in SMF-28 fiber and allowed the two fibers to be directly spliced with minimum loss. As shown in Fig. 4, OCT illumination and collection are both carried out through the single-mode (SM) core of the DCF. Meanwhile, the therapy laser is injected into the multimode (MM) end of the DCFC. The use of a single fiber greatly simplifies a future implementation in endoscopy and guarantees co-registration of the two modalities. Furthermore, an extra 1 m of SMF-28 fiber was spliced between the OCT system and the DCFC in order to reduce the ghost artifact described by Beaudette et al. [41] induced by the DCF. For the imaging head, in order to avoid damage due to the high power of the therapy laser and to maintain an achromatic behavior across the therapy and imaging lasers wavelength range, a 7 mm focal length reflective collimator (Thorlabs RC02APC-F01) was used with an UV-enhanced aluminum mirror. The collimated beam had a diameter of 2 mm before being reflected by the galvanometer (galvo) silver-coated mirror (Thorlabs GVS001) and going through the 20 mm focal length uncoated imaging lens resulting in a theoretical beam waist of 20.8 µm and a working distance of 1.06 mm in air.

Fig. 4. Setup used to obtain co-localized imaging and therapy.

Download Full Size | PDF

The therapy laser wavelength was selected because it corresponds with a peak in the absorption spectrum of hemoglobin and because it has been validated extensively in clinical applications [6,14]. Note that there is no limitation on the wavelength used for the therapy except for the DCFC, which supports wavelengths up to 1750 nm. A longer wavelength infrared laser is a more generally attractive therapy laser which would target water absorption rather than hemoglobin for a wide variety of applications. During therapy guidance, the targeted pathology will generally occupy only a portion of the OCT lateral field of view in order to provide a reference for the boundary between diseased and healthy tissue. In order to be able to selectively treat a portion of the imaging field with the therapy laser, we designed a shutter system. This was achieved by synchronizing a second galvo —the therapy galvo— to alternate between injecting the therapy laser in and out of the DCFC, as illustrated in Fig. 4. When the OCT imaging galvo was over the targeted area, the therapy galvo injected the therapy laser into the DCFC. Once the imaging galvo was no longer over the targeted area, the therapy laser was decoupled from the DCFC and sent into a beam dump by the therapy galvo. This effectively created long pulses of 10 mJ and 3 ms at 100 Hz repetition rate which corresponds to 30% of the lateral field of view. Because this setup for targeting is on the proximal end of the DCFC, it is fully compatible with endoscopic probes. A shutter was used to stop the therapy laser once the targeted therapy duration or coagulation depth was achieved. This shutter was used, instead of the therapy galvo, because the two galvos were synchronized in our setup and could not be driven independently.

In order to test our motion correction, we implemented a setup to emulate patient motion for validation. An electric motor with an unbalanced weight on its rotation axis was used to produce vibrations (see Fig. 5). A plastic airbag was used to transfer those vibrations to the sample placed on it. The airbag was wide enough to accommodate 25 cm by 15 cm samples. The frequency of the vibration could be adjusted by changing the motor frequency, and the amplitude by changing the mass of the attached weight. The vibrations reached an amplitude of around 0.3 mm at a rate able to fully decorrelate the speckle at $g^{(2)}(\tau = 10 \textrm { ms})$.

Fig. 5. Setup used to hold and vibrate sample during experiments.

Download Full Size | PDF

3.2 Software

A custom LabVIEW (National Instruments) application was developed to control both the commercial OCT system, the imaging galvo, the therapy laser galvo synchronization, and the laser shutter. The software permits three distinct modes summarized in Table 1 and described below.

Table 1. Characteristics of the three system modes

View Table | View all tables in this article

In order to study the impact of coagulation, noise, and motion on $g^{(2)}$, as well as to validate our model for motion correction, the correlation for very small values of $\tau$ was calculated. This was accomplished with the first mode called Fast mode. The imaging galvo was set to run at a 1 kHz frequency and a small scanning range while the therapy laser was left to fire continuously. The ensemble size used for calculating $g^{(2)}$ were $\Delta z=10$, $\Delta x =16$, $\Delta t=100$ in pixels or $\Delta z=77.9~\mu m$, $\Delta x=237~\mu m$, $\Delta t=100~ms$ in physical units.

In order to find the appropriate parameters for real-time monitoring and to study the validity of our approach, a second mode, called the Retrospective Diagnostic mode, was developed. It was used to develop our algorithms by recording raw data, allowing us to test our approach with different $g^{(2)}$ parameters. The Retrospective Diagnostic mode showcases how the $g^{(2)}$ can be used to add visual feedback from the therapy by providing real-time structural imaging and recording of the data for retrospective processing with the added control over the therapy laser. A shutter system allowed us to control the therapy laser exposure on the sample. Post-processing calculation of $g^{(2)}$ was performed with an ensemble size of ($\Delta z=3$, $\Delta x =10$, $\Delta t=100$) in pixels or ($\Delta z=32.5~\mu m$, $\Delta x=46.9~\mu m$, $\Delta t=100~ms$) in physical units.

The last mode, called Depth Targeting mode, consisted in the real-time monitoring of the therapy. A region of the tissue not affected by the therapy laser is manually selected to serve as a reference to quantify the bulk motion of the tissue and apply motion correction. An ensemble of pixels is manually selected at the targeted depth and used to check if the depth has been reached in that region. By calculating $g^{(2)}$ in selected regions, computation in real time becomes possible. For a field of view of 2.34 mm (300 pixels) in depth and of 25 mm (800 pixels) laterally, the ensemble size used was ($\Delta z=3$, $\Delta x=50$, $\Delta t=10$) in pixels or ($\Delta z=32.5~\mu m$, $\Delta x=234.4~\mu m$, $\Delta t=100~ms$) in physical units in air. The ensemble size was chosen to have a high resolution in $z$ and $t$ while keeping a good estimation of $g^{(2)}$ in order to have more accurate automatic coagulation depth determination. It was observed that this larger ensemble size allowed for a more accurate motion correction. Since our goal was coagulation depth with no need to control the lateral extent of the coagulation, a larger size in $x$ was used to allow us to have a finer resolution in $z$.

To perform real-time motion correction, two reference ensembles on both sides of the treated area were used to calculate $g^{(2)}_{nc}$, using the method described in Sec. 2.3. The two ensembles are concatenated into one to calculate a $g^{(2)}_{nc}ref$ value over a larger ensemble, and improve the accuracy of the motion estimation. The ensemble size were $\Delta z =150$, $\Delta x =100$, $\Delta t=10$. We established whether the targeted coagulation depth was reached by testing for two conditions:

(14)$$\begin{aligned}& C_1:{ \textbf{if} \quad} g^{(2)}_{nc\&mc}(\tilde{t}) < \xi_1 \textrm{ for any }\tilde{t} < t. \\ & C_2:{ \textbf{if} \quad} g^{(2)}_{nc\&mc}(t) > \xi_2 {\textbf{ and if} }C_1 \textrm{ is satisfied}, \end{aligned}$$

where $t$ is the current time index at which the conditions are checked, $\tilde {t}$ is any time index before the current one and $\xi _1$ and $\xi _2$ are two threshold values. $C_1$ consists of checking if the selected region has started to coagulate. As the tissue begins to coagulate, $g^{(2)}$ drops. We can therefore set a certain threshold ($C_1$) indicative of the start of the process. As coagulation nears completion, $g^{(2)}$ will then return to 2 as the tissue becomes static/stable again. $C_2$ is set as the threshold indicating the completion of the coagulation process. The values of $\xi _1$ and $\xi _2$ were both determined experimentally, as explained in Sec. 4.2. The system was capable of verifying both conditions at a faster rate than the B-scan frequency. The data process in summarized in Fig. 6. The SNR was estimated using the $R$ versus signal intensity curve shown in Fig. 7. When $\tilde {g}^{(2)}_{nc} < 1.3$ (that is, motion produced a strong decorrelation of the signal), the conditions $\{C_1,C_2\}$ are not computed until $\tilde {g}^{(2)}_{nc} \ge 1.3$. This allowed us to perform coagulation depth tracking despite sudden large motion.

Fig. 6. Data process steps for real-time monitoring.

Download Full Size | PDF

Fig. 7. Estimated SNR vs. noise intensity calibration curve.

Download Full Size | PDF

3.3 Histology

The rat tongue and abdominal muscle samples were taken from a rat a few hours post mortem. The samples were enveloped in a hydrated gauze with Hank’s Balanced Salt Solution (HBSS) during transport and kept on an hydrated gauze during the experiment, as shown by the Fig. 5. We didn’t soak the sample in HBSS in order to avoid the blood in the tissue to diffuse in the solution. Since we are targeting hemoglobin absorption, it is essential to preserve the blood content of the sample.

Once the experiment was finished, the optimum cutting temperature compound was used for frozen section histological processing. The sample was then stored on dry ice until being brought for processing. Hematoxylin and eosin (HE) histology with nitroblue tetrazolium chloride (NBTC) stains for the activity of nicotinamide adenine dinucleotide diaphorase (NADH-diaphorase) were used on slices 10 µm thick in a fashion similar to [29,57].

4. Results

4.1 Noise and motion correction

In a first step, we estimated the SNR calibration curve by calculating $g_{nc}^{(2)}$ in 44 B-scans acquired while imaging rat muscle tissue before the start of the therapy. The expected, relatively linear relationship between $R_j$ and $I_j$ is presented in Fig. 7. Given Eq. (10) predicts a linear proportionality between the SNR and $I$, a simpler calibration which simply finds $\langle N \rangle$ could work equally well.

The estimation fails at very low SNR values as most of the input is noise. As for other SNR levels, the difference can be explained as well. First, the limited size of ensembles, which adds some uncertainty to the calculated values. Second, the assumption that the noise floor is constant and independent of $z$, $x$ or $t$. Third, in a hypothetical incoherent cross-sectional image —that is in the absence of speckle— the sample reflectivity will not be completely uniform inside the area corresponding to a single ensemble due to the heterogeneous tissue structure. This implies that assigning a single average intensity to an ensemble is an approximation. Finally, estimating the SNR on very bright spots is limited by saturation of the signal. Despite these limitations, the SNR estimation follows the predicted shape in Eq. (12) and allowed us to correct for noise as shown in Fig. 8.

Fig. 8. $g^{(2)}$ before (circle) and after (asterisk) noise correction for low SNR, $R = 9.50 dB$, (dark blue) and high SNR, $R = 15.91 dB$, (light blue) values. Color palette from [58].

Download Full Size | PDF

The next step was to measure $g_{nc}^{(2)}$ as a function of $\tau$, with the Fast mode described in Sec. 3.2, to determine the optimal value for monitoring the coagulation process. Equation (11) was used with a fixed value of $z$, $x$ and $t$. This corresponds to using the autocorrelation function to study a static process, because we will limit the analysis to time intervals in which the sample can be assumed to be in a single state (static, coagulating or after coagulation).

In order to observe the effect of the noise correction, the coagulation was performed on an excised sample of rat muscle containing enough hemoglobin to absorb light from the therapy laser. Imaging was performed at a 1 kHz B-scan rate, providing us with a resolution of 1 ms for $\tau$ —the shortest decorrelation measurable— and a resolution of $100$ ms for $t$ considering the ensemble size $\Delta t$ —the time scale over which the tissue is assumed to be in the same state— over a width of 0.35 mm, while the therapy laser was left continuously on with an average power of 1.0 W at the sample. The $g^{(2)}$ was calculated at two different depths ($z$ values), for different SNR levels, $R = 8.91$ and $R = 39.03$, with a fixed $t$ for each $z$ during which the tissue was undergoing coagulation at that depth. Figure 8 shows how noise correction impacts the calculated $g^{(2)}$ for two different SNR values. As one would expect, the correction has a more significant impact for lower SNR values. The difference is also more pronounced at shorter delays $\tau$ (or higher $g^{(2)}$ values). Implementing noise correction will thus allow for speckle correlation measurements deeper within the tissue, as the OCT SNR dramatically diminishes with depth and coagulation state. This is important for our Depth Targeting mode to work properly. The drop of $g^{(2)}$ between $\tau = 0$ and $\tau = 1$ ms is well-known and explained by the noise pattern becoming decorrelated between B-scans, which is not the case for $\tau = 0$, since we are correlating the same noise pattern.

In our next step, we first validated our simulation model by comparing the total correlation for a range of coagulation-induced and motion-induced decorrelation values for both simulated and experimental data. The process to create the simulated data is described in Sec. 7. For the experimental data, we measured the coagulation of rat muscle at 1 kHz with 1.0 W therapy laser power, with a $\tau$ resolution of 1 ms. Figure 9(a) corresponds to this comparison; the horizontal axis shows the value of $g^{(2)}$ for a range of coagulation-induced decorrelation, while the vertical axis shows the combined effect of coagulation- and motion-induced decorrelation. The dashed lines with unity slope represent simulated and experimental $g^{(2)}$ with no motion, which overlap in Fig. 9(a). Motion-induced decorrelation was added by digitally shifting the second correlated ensemble by 1 pixel (7.8 µm) to both the simulated and experimental data (solid lines). The $g^{(2)}$ value drops very similarly in both cases. This validated our approach to create the lookup table for motion correction, which was based on simulated data with digitally added motion. We note that the simulated data was noise free and that the experimental data was corrected for noise.

Fig. 9. (a) Validation of the simulation and motion emulation. $g^{(2)}$ of speckle for simulated (yellow) and experimental (orange) tissue coagulation. Dashed lines represent the data without added motion, both lines are overlapping. Solid lines show $g^{(2)}$ with motion-induced decorrelation added to the coagulation-induced decorrelation. Motion was induced by digitally shifting one of the ensembles by 1 pixel. The figure confirms that motion-induced decorrelation affect simulated and experimental coagulation-induced decorrelation in the same way. (b) Experimental results of the motion correction procedure. $g^{(2)}$ during coagulation (pink curves) and after full coagulation (purple curves), with (dashed curves) and without (solid curves) motion correction. Color palette from [58].

Download Full Size | PDF

Next, the effect of motion correction on experimental data was evaluated. The OCT scan rate was set at 1 kHz over a small field of view and the therapy laser power was set to 1.0 W. Motion was induced with a motor, as described in Sec. 3 $g^{(2)}$ was computed for a fixed $z$ and $x$ for two fixed $t$ corresponding to during and after coagulation with the laser still on in both instances. Coagulation was assessed visually after laser exposure.

In this fast mode, there is no untreated region of the sample that can be used for determining $\tilde {g}^{(2)}_{nc}$ during the treatment. Instead, we assumed that the magnitude of the motion would remain approximately constant during the entire experiment. We thus estimated the motion by determining $\tilde {g}^{(2)}_{nc}$ before the therapy started, and used this same value for all subsequent motion correction. The motion was estimated to range from 6.4 to 6.8 µm per 10 ms. Results are shown in Fig. 9(b).

Similarly to the noise correction, motion correction makes the most difference for speckle that should be correlated without motion. Thus, the correction allows us to confirm whether a speckle pattern is stable in the presence of motion. Figure 9 also shows the correlation of the speckle pattern during and after coagulation. As expected, the $g^{(2)}(\tau )$ value decreases much faster during coagulation than after.

4.2 Distinguishing coagulated from coagulating tissues

As we see with Fig. 9(b), coagulated tissue can be distinguished from tissue still undergoing coagulation. In this section, we study the variation in time of $g^{(2)}$, for a different tau values. This allows the identification of the optimal tau value, for which the contrast between coagulating and fully-coagulated tissue is maximized. With this in mind, we monitored the coagulation of abdominal rat muscle exposed to 10 mJ pulses at 100 Hz —which is the nominal setting in our system for simultaneous targeted therapy and monitoring— during 3 s. With our system in Retrospective Diagnostic mode (described in Sec. 3.2), the therapy laser was put in pulsed mode as to fire only at the center of the OCT field of view with a pulse energy of 10 mJ at 100 Hz, i.e. using 1.0 W of average power.

Figure 10 shows the dynamic $g^{(2)}$ as a function of time $t$ during the exposure of the sample to the therapy laser for four different values of $\tau$. $g^{(2)}$ was calculated for an area at the center of the field of view, close to the surface, where the SNR is the highest and the coagulation was guaranteed to complete. In this case, a pixel size ensemble of $z=3$, $x=50$, and $t=10$ was used for a more accurate estimation of $g^{(2)}$ and since high spatial resolution was not a concern. The coagulation spot was visible at the end of the therapy, confirming that tissue did indeed coagulate. From left to right, vertical dashed lines represent moments when the laser was turned on, when tissue motion due to thermal expansion stops (evaluated visually), and when coagulation is reached. For all $\tau$ values, a drop in the $g^{(2)}$ value is observed shortly after the laser is turned on. There is often a sudden thermal expansion of the tissue during the initial phase of the coagulation. The expansion is easily identifiable visually. We can see the end of this expansion at the second line when $g^{(2)}$ increases slightly for all values. Finally, at the third line, we can see an increase in $g^{(2)}$ for all values of $\tau$. We estimate the tissue to have finished coagulating at this moment. These results were used in determining our choice of threshold values and $\tau$ for the real-time monitoring, which were set from inspection of Fig. 10 to $\tau =10$ ms, $\xi _1=1.4$ and $\xi _2=1.9$.

Fig. 10. $g^{(2)}$ for $\tau$ ranging from 10 ms to 40 ms (from dark to pale blue) during laser exposure. From left to right, vertical dashed lines indicate when the laser is turned on, the end of thermal expansion, and when we estimate coagulation to be complete. Color palette from [58]

Download Full Size | PDF

The use of threshold values is, however, not necessary to visualize the progression of the thermal coagulation area of the tissue. Figure 11 shows a superposition of the decorrelation map over the OCT image for the same experiment. Ensemble size was returned to ($z=3$, $x=10$, $t=10$), since spatial resolution was necessary to know which part of the tissue was affected by coagulation. From left to right, Fig. 11 shows the OCT intensity tomogram (a), the $g^{(2)}$ map without (b) and with (c) noise correction and the $g^{(2)}$ overlay in color over the OCT tomogram (d). The difference made by noise correction can be seen by comparing the $g^{(2)}$ map without (b) and with (c) noise correction. The red bars shows the depth at which the system detects undergoing coagulation, $g^{(2)} < \xi _1$, due to noise in the untreated region of the tissue. The noise correction thus extends the monitoring range by approximately 0.5 mm. The overlay preserves the structural information from OCT while highlighting the progression of the coagulation. Since Labview does not permit creating the overlay, it was performed in retrospective processing. The refractive index of rat muscle was assumed to be $n=1.39$ for adjusting the scale of the $z$ axis [59,60]. Coagulation can be seen as the transition from correlated speckle to decorrelated and back to correlated. Coagulated tissue is thus visible as a correlated area in yellow over a decorrelated arc in blue. The progression of the thermal coagulation can also be observed in Fig. 12 and Visualization 1. Coagulation overlay was not calculated with Eq. (14) since we found the threshold method to not be reliable with limited ensemble sizes.

Fig. 11. Processing steps during Retrospective Diagnostic mode. (a) OCT structural image during therapy, (b) corresponding $g^{(2)}_{n}$, (c) corresponding $g^{(2)}_{nc}$, and (d) overlay of the OCT intensity (encoded in the luminance) and $g^{(2)}_{nc}$ encoded in the hue of a perceptually uniform isoluminant color map [61]. The red bars in (b) and (c) indicate the depth at which the system would have incorrectly detect undergoing coagulation, $g^{(2)} < \xi _1$, due to noise in the untreated region. Scale bars are 0.5 mm.

Download Full Size | PDF

Fig. 12. Therapy monitoring overlay highlighting the progression of the coagulation using the correlation of the speckle pattern after 0, 0.5, 1.0, 2.0, 3.0 and 4.0 seconds of exposure by 10 mJ pulses therapy laser at 100 Hz. Scale bars are 0.5 mm.

Download Full Size | PDF

Figure 12 shows a sequence of OCT images during the therapy of the sample at different time points. The periodic horizontal lines in the intensity are a product of birefringence of muscle at baseline —which is affected by the coagulation process— and the absence of a polarization-diverse receiver in the OCT system. Homogeneous tissue birefringence produces a regular change in the polarization state of light with depth. Given that the single reference arm polarization is fixed, interference between sample and reference arms follows periodic oscillations as a function of depth, essentially reproducing Malus’s Law [62]. Dark bands correspond to depths at which tissue birefringence causes backscattered light to have perpendicular polarization with respect to that of light from the reference arm.

The final experiment in Retrospective Diagnostic mode consisted of the coagulation monitoring of a moving sample and confirming our result though histology. The parameters used were the same, 10 mJ pulses at 100 Hz, however duration was 5 seconds instead, to adjust for the different optical and thermal properties of the sample. B-scan rate was the same 100 Hz and a $\tau =10$ ms was used to calculate $g^{(2)}$. Motion was introduced using a vibrating motor and quantified using ensembles on both sides of the targeted area as explained in Sec. 3.2. Figure 13 shows the coagulation monitoring of the tongue of a rat subjected to motion before (a) and after (b) motion correction, as well as the corresponding histology section (c). The OCT intensity image is shown in (d), which does not exhibit any clear indication of the coagulation depth. Images shown were acquired post laser exposure since we observed that the coagulation continues for a few tens of milliseconds after the laser is turned off. Despite the expected change in shape of the tissue after histological processing —particularly the fold of the tongue which became more pronounced— the coagulation depth observed in OCT matches fairly well that measured with histology. Motion correction improves the visibility of the correlated area corresponding to fully coagulated tissue. More importantly, noise- and motion-correction together are essential in the use of threshold values —described in Eq. (14)— for the real-time monitoring of the coagulation depth and automatic interruption of the therapy laser presented in the next section.

Fig. 13. Coagulation monitoring overlay in rat tongue in presence of vibration (a) without motion correction and (b) with motion correction compared to (c) the histology section and (d) the OCT structural image. Red markers of identical lengths (0.55 mm) are placed on each image to facilitate comparison. Scale bars are 0.5 mm.

Download Full Size | PDF

4.3 Monitoring in real-time

The previous images were all created in post-processing with the Retrospective Diagnostic mode. Although images shown in this section were created in post-processing as well, the Depth Targeting mode was used for which a desired coagulation depth is defined interactively on the initial structural image. The correlation is calculated at this targeted depth in real time to perform the coagulation-threshold values checks, outlined Eq. (14).

We were able to detect the moment the coagulation reached the desired depth and shut the laser off. The experimental parameters are detailed in Table 1. Figure 14 shows the results from two real-time experiments which were performed targeting two different depths. The overlays correspond to the state of the tissue a few tens of milliseconds after the laser was turned off since the coagulation process continues after exposure. We can see that in the case of the shallower targeted depth, we seemed to have coagulated deeper than intended. This might be due to the coagulation continuing after the therapy laser closed. In a future implementation, the increase in coagulation depth after the laser is off can be determined experimentally, and the laser can be shut off before in order to reach the targeted depth once coagulation ceases. The difference was not large for the deeper target depth since we observed that the post laser shut off coagulation slows down deeper in tissue.

Fig. 14. Experimental results in which laser therapy depth coagulation was automatically stopped in real time (a–b) at 0.48 mm and (c–d) 1.04 mm below the tissue surface, indicated by the red line. The therapy monitoring overlay highlighting coagulation is shown (b–d) on the right while the corresponding OCT structural image is shown (a-c) on the left. Scale bars are 0.5 mm.

Download Full Size | PDF

We performed a final experiment in the presence of induced motion using the technique described previously, with the corresponding results in Fig. 15. We targeted a specific depth and compared the real-time estimation of coagulation depth from our technique with histology. The overlays show the tissue tens of milliseconds after the laser was turned off. The OCT overlays shown correspond to only noise-corrected decorrelation and noise- and motion-corrected. The targeted depth, marked by a red line, corresponds well to the coagulation depth determined in the histology. The full video of the results of this experiment are presented in Visualization 2.

Fig. 15. Therapy monitoring overlay highlighting coagulation during laser therapy (a) without motion correction and (b) with motion correction with (c) histology and (d) OCT structural image for comparison. The therapy was automatically stopped in real time when the targeted depth of 0.94 mm under the surface was reached, as shown by the red line. Scale bars are 0.5 mm.

Download Full Size | PDF

We quantified the motion by obtaining the vertical displacement required to get the drop in correlation observed in $\tilde {g}^{(2)}_{nc}$ from the lookup table. We estimate that the decorrelation induced by the motion corresponds to a 55 µm displacement in the axial direction or 108 µm in the lateral direction over 10 ms.

5. Discussion

With our laser therapy and OCT system, we demonstrated how the normalized second-order autocorrelation function $g_n^{(2)}$ of the OCT intensity signal could be used to monitor thermal coagulation depth. Noise correction on $g^{(2)}$ gives our system an increased monitoring depth while motion correction makes it more robust to bulk tissue motion. Furthermore, we showed that it was possible to verify whether the coagulation had reached the targeted depth in real time. Our results were confirmed on an ex-vivo sample (rat muscle) through NBTC histology. Since monitoring is based on the signal intensity correlation, it does not require the phase signal to be stable. In particular, the compute-intensive phase stabilization and complex-amplitude tomogram interpolation required to correct for non-uniform rotation distortion in endoscopic systems [30] is not required, and instead a simpler intensity-based approach would be sufficient [63]. Combined with the use of a single DCF coupled to a DCFC, the implementation of this approach into clinical catheter systems is greatly simplified and the co-registration of the therapy and imaging are guaranteed.

The coagulation progression was monitored in two ways: one aimed at visualizing the tissue region under active coagulation and letting the operator interpret the corresponding treatment depth (Retrospective Diagnostic mode); the other aimed at automatically determining the treatment depth for automatic laser shut down once the target depth is reached (Depth Targeting mode). In both cases, the depth was monitored over the span of several hundred microns corresponding to the typical thickness of the affected tissue layers in potential applications of our technique: 200–500 µm for the epithelium in the esophagus, 1 mm for respiratory papillomatosis and 0.35 mm for diabetic retinopathy.

The Retrospective Diagnostic mode was implemented in retrospective processing for our experiments, however it was designed for and is fully compatible with a real-time implementation. Visual feedback is provided by means of overlaying the hue-encoded speckle correlation value and the luminance-encoded OCT structural image, preserving the structural information while adding coagulation contrast during the treatment. The progression of the coagulation is readily seen by the transition from a decorrelated state, when tissue is being heated, to correlated, once it has reached coagulation. In this case, we did not automatically detect the coagulation as proposed in Eq. (14). The reason for this was the small ensemble size used, ($z=3$, $x=10$, $t=10$) pixels. We found the threshold method to not to be reliable with limited ensembles. In this case, visual estimation by the operator is straightforward. A neural-network approach could potentially replace the operator making the coagulation depth monitoring more robust in this mode.

For the moment, the Retrospective Diagnostic mode was implemented in post-processing due to LabVIEW not allowing the creation of an overlayed image and its limited $g^{(2)}$ computing rate. We confirmed that our custom LabVIEW program using only CPU computing (no GPU) was capable of calculating $g^{(2)}$ on a 400x800 pixels frame using ensembles comprising 10 B-scans at a rate of 5 Hz with both noise and motion corrections. We expect that the use of faster programming language and the implementation of GPU-computing would easily allow our method to run in real-time.

The depth targeting mode was performed in real time, automatically shutting off the therapy laser when the coagulation had reached the targeted depth. We reduced the computation and display demands by calculating the correlation for a single ensemble situated at the target. To increase the reliability of the shutting off of the laser, we used a larger ensemble ($z=3$, $x=50$, $t=10$).

Despite high success, in some rare cases a sudden thermal expansion of tissue caused a false positive in the depth detection. The thermal expansion resulted in motion of the untreated, deeper part of the tissue, causing a drop in $g^{(2)}$ followed by an increase —after expansion stops— and thus detected as complete coagulation at a depth not yet reached by the treatment. This phenomenon, easily identifiable visually, caused the automatic stop to trigger, ending the therapy prematurely could be automatically detected by calculating the cross-correlation of the OCT signal in depth to detect large-scale axial motion of the tissue and disable momentarily the coagulation detection. Another limitation of the depth targeting mode is that there is a minimum depth at which we can stop the coagulation, roughly 100 µm in rat muscle tissue. This is due to the laser power spread in depth which causes a coagulation front of finite depth.

In terms of the calculation of the correlation, we observed that the initial noise calibration of the system remained accurate despite changes in the tissue scattering properties, implying that the noise floor of the system is indeed constant throughout the therapy. Thanks to the noise correction, the monitoring depth was limited by the power of our therapy laser and not by the monitoring technique. With a different therapy approach able to deposit energy deeper in tissue the monitoring depth could become limited by the OCT signal attenuation.

Regarding our motion correction, we found it was essential for the depth targeting mode. An obvious limitation of the method arises when the motion causes the speckle pattern to become fully decorrelated, as there is no correction possible to recover the value of the correlation in absence of motion. To mitigate this, we limited the correction of motion-induced decorrelation to $g^{(2)} > 1.3$; for $g^{(2)} \leq 1.3$ we consider any measure of $g^{(2)}$ to be inaccurate and the checking of the coagulation conditions of Eq. (14) are not performed. This allows the mode to continue working despite sudden and large motion. Furthermore, a warning message is displayed when this occurs. For our system configuration, this corresponded to approximately 109 µm of axial displacement or 246 µm of lateral displacement between two B-scans. It is interesting to note that the dependence of speckle size on the system resolution implies that a system with poorer resolution would be capable of correcting for greater displacement, at the expense of the spatial resolution and the contrast between the coagulation front and coagulated tissue.

An important aspect of our approach was the determination of the optimal value of $\tau$ to maximize the contrast between the coagulated area and static tissue, as well as the optimal threshold values $\xi _1$ and $\xi _2$ for the automatic depth targeting. These values are expected to change depending on the absorption of therapy light, tissue thermal properties and the therapy laser power, wavelength and mode —CW or pulsed [64]. Therefore, prior testing is required to find the optimal parameters for each application.

Another important consideration is the optimal size for the ensemble used to calculate $g^{(2)}$. As stated in the theory, larger ensembles have the advantage of estimating the correlation with a greater reliability, at the cost of a lower resolution in space and/or time. Furthermore, we found that motion correction necessitates larger ensembles to properly estimate and correct the motion. Our choice of ensemble size provided us with a resolution of $\sim 32.5$ µm in $z$ and 100 ms in time for the automatic coagulation depth targeting mode.

Improvements to the experimental setup include custom coating to the imaging lens to improve transmission for both the OCT and therapy laser transmission, especially in the context of endoscopic imaging. Our proof of principle relied on a chopped CW laser, which greatly limited the available power. In a next iteration, we envision using a pulsed laser in a wavelength range targeting an absorption peak of water.

Furthermore, the fact that the random movement of scatterers mimic the shape of the autocorrelation function of coagulating tissue opens an alternative avenue for correcting for motion. It is known that diffusive motion, another example of scatterer random motion, is described by an exponential decay of the autocorrelation function, in contrast to the Gaussian shape of bulk motion [65–67]. By calculating $g^{(2)}$ as a function of $\tau$ for a given $t$, it is possible to fit an exponential and a Gaussian function, in similarity with the approach to decouple diffusion- from translation-induced decorrelation in dynamic light scattering, and distinguish coagulation- from motion-induced decorrelation. Preliminary testing [68] shows great potential, but requires more computational power to calculate the autocorrelation function as a function of $\tau$ in real time.

Future work translating our technique to in-vivo experiments would require addressing some remaining challenges. Living tissue would introduce further sources of decorrelation besides coagulation and motion, such as blood flow inside blood vessels, potential bleeding and edema. However, decorrelation effects from these are expected to manifest on time scales significantly slower than coagulation-induced decorrelation. Their impact on our method would be small; optimal threshold values $\xi _1$ and $\xi _2$ might need to be optimized. Motion correction might present new challenges as well. Our correction method expects uniform bulk motion across B-scans and does not account for heterogeneous distortions. A potential solution is the use of tracking of tissue features (as opposed to speckle) using first, for instance, the speeded up robust feature (SURF) algorithm [69], followed by the Kanade-Lucas-Tomasi (KLT) feature tracking algorithm [70] across each pair of B-scans. This would produce a map of large-scale tissue deformation that would be used to suppress motion-induced decorrelation. There are real-time implementations of these algorithms [71]. Finally, further studies in vivo would require the use of endoscopic probes: an advantage of our use of DCFCs is that replacing a standard single-mode fiber with a DCF is relatively straightforward and would be enough for implementing this technique endoscopically, as long as the therapy beam spot size is of adequate size for coagulation [41].

6. Conclusion

In this work we presented a new laser therapy monitoring technique that allowed us to monitor the progression of thermal coagulation of excised animal tissue. A 532 nm therapy laser was used with OCT imaging through a single double-clad fiber compatible with endoscopic use. The use of a single fiber allows for simultaneous and co-registered imaging and therapy and greatly facilitates its implementation in endoscopic probes. We present the use of a normalized second-order autocorrelation function to measure the decorrelation of the speckle pattern occurring during coagulation. We implement an analytical noise correction that extends the therapy monitoring depth beyond the penetration depth of the therapy laser. We also implement a numerical motion correction technique that makes our technique robust against bulk tissue motion. Using a threshold scheme, we are able to detect, in real-time, when the coagulation had reached the targeted depth and automatically turn the therapy laser off. Relying on the dynamics of speckle during coagulation addresses fundamental limitations of techniques that rely on the change of tissue birefringence —which require tissue with birefringence at baseline and a polarization-sensitive OCT system, or the attenuation of the OCT signal —which are known to underestimate the therapy depth. Finally, in contrast with alternative methods that rely on the complex-amplitude signal correlation, our monitoring scheme does not require the phase of the OCT signal, further facilitating real-time implementation in endoscopic OCT systems, thanks to the enabling of simplified correction approaches for non-uniform rotation distortion artifacts.

7. Appendix A - Speckle simulation

In this Appendix we present the procedure used to numerically simulate speckle patterns to study the combined effect of motion and thermal coagulation on $g^{(2)}$. We randomly distributed an ensemble of scatterers and simulated light from an illumination beam being reflected off of those scatterers, producing a two-dimensional speckle pattern [56]. Coagulation was simulated by randomly moving the scatterers individually, while bulk tissue motion was added with a sub-pixel displacement. $g^{(2)}$ calculated from the resulting speckle patterns were assembled into the lookup table presented in Fig. 3. Simulation parameters are included in Table 2.

The first step consisted of simulating the position of randomly distributed and moving scatterers. We start by defining our simulated environment pixels $z$ and $x$ positions, $p_{z0}$ and $p_{x0}$, as follows:

(15)$$\begin{aligned}p_{z0}(z) = (z-1) \cdot px \textrm{ for } z \in \{0,1,2,\dots,s_z-1\}, \\ p_{x0}(x) = (x-1) \cdot px \textrm{ for } x \in \{0,1,2,\dots,s_x-1\},\end{aligned}$$

where $z$ and $x$ are the pixel indexes, $s_z$ and $s_x$ are the size of the simulated environment and $px$ is the physical size of the pixels, assumed to be the same in $z$ and $x$ for simplicity. We then randomly distributed the scatterers:

(16)$$\begin{aligned}p_{zs}(z,x,n,t=0) = \epsilon_{1z}(z,x,n) \cdot px + p_{z0}(z) , \\ p_{xs}(z,x,n,t=0) = \epsilon_{1x}(z,x,n) \cdot px + p_{x0}(x) ,\end{aligned}$$

where $p_{zs}(z,x,n,t=0)$ and $p_{xs}(z,x,n,t=0)$ are the physical locations, in the $z$ and $x$ directions respectively, of the $n$-th scatterer in pixel ($z$, $x$). There is thus multiple simulated scatterers per pixel for a total number of $s_n$ scatterers per pixel. $\epsilon _{1z}$ and $\epsilon _{1x}$ are randomly generated numbers uniformly distributed from 0 to 1. Random motion is then added to simulated tissue coagulation [68]

(17)$$\begin{aligned}p_{zs}(z,x,n,t) = p_{zs}(z,x,n,t-1) + m_a \cdot \epsilon_{2z}(z,x,n,t), \\ p_{xs}(z,x,n,t) = p_{xs}(z,x,n,t-1) + m_a \cdot \epsilon_{2x}(z,x,n,t),\end{aligned}$$

where $\epsilon _{2z}(z,x,n,t)$ and $\epsilon _{2x}(z,x,n,t)$ also are randomly generated number with a uniform distribution between 0 and 1. The amplitude of the movement is scaled by changing $m_a$ value. We then simulate how OCT light is reflected by our scatterers. We simulated a Gaussian point spread function (PSF) centered on each pixel and summed the resultant electric field reflectivities

(18)$$\begin{aligned}E(z,x,t) = \sum_{Z=1}^{s_z} \sum_{X=1}^{s_x} \sum_{N=1}^{s_n} \exp{\left( \frac{-[z-p_{zs}(Z,X,N,t)]^2-[x-p_{xs}(Z,X,N,t)]^2}{2 w^2} \right)} \\ \exp{\left({-}2 i \pi p_{zs}(Z,X,N,t) \right)},\end{aligned}$$

where $E$ is the resulting complex electric field on each pixel and $w$ is the beam waist of the PSF assuming Gaussian illumination. The simulation thus accounted for both the amplitude and phase of the backscattered light. We are, however, interested in the signal intensity:

(19) $$I(z,x,t) = |E|^2.$$

The process in Eqs. (17) and (18) is repeated to obtain a speckle pattern changing in time in response to the dynamics of the scatterers. The next step consisted of calculating $g^{(2)}$. We defined the ensembles $I_1$ and $I_2$ as follows:

(20)$$\begin{aligned}I_1 &\equiv \{I(z,x,t=0)\}, \forall z \textrm{ and } \forall x \textrm{ and} \\ I_2(\tau) &\equiv \{I(z,x,t=\tau)\}, \forall z \textrm{ and } \forall x. \end{aligned}$$

The variable $\tau$ defines the delay. Since we also wanted to simulate bulk tissue motion, instead ot adding it in Eq. (17), we decided to shift $I_2(\tau )$ by a subpixel amount based on previous works on subpixel registration [72–74]:

(21)$$\begin{aligned}&Z = \mathcal{F}(I_2(\tau)), \\ &z_{shift} = \exp{({-}2i \pi \Delta)}, \\ &Z_{shifted} = Z \cdot z_{shift}, \\ &I_2(\tau,\Delta) = \mathcal{F}^{{-}1}(Z_{shifted}), \end{aligned}$$

where $\mathcal {F}$ denotes a two-dimensional discrete Fourier transform implemented through a fast Fourier transform algorithm, and $\mathcal {F}^{-1}$ is the corresponding inverse discrete Fourier transform. $\Delta$ is the shift value in pixels ranging from 0 to $\Delta _{\textrm {max}}$. Since we shifted $I_2(\tau ,\Delta )$ towards the bottom, the top of the speckle pattern will no longer be spatially correlated. We thus cropped the wrapped regions before calculating g2. $I_1$ and $I_2(\tau ,\Delta )$ become

(22)$$\begin{aligned}I_{1c} \equiv \{I_1(z,x)\} \ \forall z \in \{\Delta_{\textrm{max}}+1,\Delta_{\textrm{max}}+2,\dots,s_z\} \textrm{ and } \forall x \textrm{ and} \\ I_{2c} \equiv \{I_2(\tau,\Delta)(z,x)\} \ \forall z \in \{\Delta_{\textrm{max}}+1,\Delta_{\textrm{max}}+2,\dots,s_z\} \textrm{ and } \forall x.\end{aligned}$$

We finally calculated $g^{(2)}$ for varying values of $\tau$ and $\Delta$

(23)$$g^{(2)}(\tau,\Delta) = g^{(2)}_n (I_{1c},I_{2c}(\tau,\Delta)),$$

where $g^{(2)}_n$ is defined in Eq (5). To calculate the decorrelation map, we used the parameters in Table 2.

Table 2. Speckle simulation parameters

View Table | View all tables in this article

Funding

National Institutes of Health (K25 EB024595); Fonds de recherche du Québec – Nature et technologies (2018_PR_208631).

Acknowledgments

The authors would like to acknowledge Édith Ducharme, Xavier Attendu and Martin Poinsinet De Sivry-Houle for useful discussions, and Irene Londoño for histology.

Disclosures

NUP: Terumo Corporation (P), CB: Castor Optics, inc. (I,P).

References

1. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]

2. S. Gabay, “Highlights of laser-tissue interaction mechanism,” in Laser Florence 2000: A Window on the Laser Medicine World, vol. 4606 (International Society for Optics and Photonics, 2001), pp. 30–33.

3. W. Sheng, S. He, W. J. Seare, and A. Almutairi, “Review of the progress toward achieving heat confinement–the holy grail of photothermal therapy,” J. Biomed. Opt. 22(8), 080901 (2017). [CrossRef]

4. L. Polese, I. Angriman, M. Scarpa, D. Pagano, P. Parente, F. Erroi, M. Frego, D. F. D’Amico, and L. Norberto, “Diode laser treatment of Barrett’s esophagus: long-term results,” Lasers Med. Sci. 26(2), 223–228 (2011). [CrossRef]

5. K. K. Wang and J. Y. Kim, “Photodynamic therapy in Barrett’s esophagus,” Gastrointest. Endosc. Clin. North Am. 13(3), 483–489 (2003). [CrossRef]

6. L. Gossner, A. May, M. Stolte, G. Seitz, E. G. Hahn, and C. Ell, “KTP laser destruction of dysplasia and early cancer in columnar-lined Barrett’s esophagus,” Gastrointest. Endosc. 49(1), 8–12 (1999). [CrossRef]

7. “Early Photocoagulation for Diabetic Retinopathy,” Ophthalmology 98, 766–785 (1991).

8. “Photocoagulation for Diabetic Macular Edema: Early Treatment Diabetic Retinopathy Study Report Number 1 Early Treatment Diabetic Retinopathy Study Research Group,” Arch. Ophthalmol. 103, 1796 (1985).

9. S. Sivaprasad and G. Dorin, “Subthreshold diode laser micropulse photocoagulation for the treatment of diabetic macular edema,” Expert Rev. Med. Devices 9(2), 189–197 (2012). [CrossRef]

10. C. R. Munnerlyn, “Lasers in opthalmology: Past, present and future,” J. Mod. Opt. 50(15-17), 2351–2360 (2003). [CrossRef]

11. S. S. Hayreh, P. A. Podhajsky, and M. B. Zimmerman, “Natural History of Visual Outcome in Central Retinal Vein Occlusion,” Ophthalmology 118(1), 119–133.e2 (2011). [CrossRef]

12. “Argon Laser Scatter Photocoagulation for Prevention of Neovascularization and Vitreous Hemorrhage in Branch Vein Occlusion: A Randomized Clinical Trial,” Arch. Ophthalmol. 104, 34 (1986).

13. M.-L. Kuet and M. J. Pitman, “Photoangiolytic Laser Treatment of Recurrent Respiratory Papillomatosis: A Scaled Assessment,” J. Voice 27(1), 124–128 (2013). [CrossRef]

14. F. Benboujja, S. Bowe, C. Boudoux, and C. Hartnick, “Utility of Optical Coherence Tomography for Guiding Laser Therapy Among Patients With Recurrent Respiratory Papillomatosis,” JAMA Otolaryngol. Head & Neck Surg. 144(9), 831–837 (2018). [CrossRef]

15. J. T. Wei, P. Sharma, and N. J. Shaheen, “Endoscopic ablation therapies for Barrett’s esophagus: a meta-analysis of stricture complications,” Am. J. Gastroenterol. 98 (2003).

16. N. J. Shaheen and D. J. Frantz, “When to consider endoscopic ablation therapy for Barrett’s esophagus,” Curr. Opin. Gastroenterol. 26(4), 361–366 (2010). [CrossRef]

17. S. Eluri and N. J. Shaheen, “Endoscopic Eradication Therapy in Barrett’s Esophagus,” Tech. Gastrointest. Endosc. 19(3), 137–142 (2017). [CrossRef]

18. J. Shah, S. R. Aglyamov, K. Sokolov, T. E. Milner, and S. Y. Emelianov, “Ultrasound imaging to monitor photothermal therapy – Feasibility study,” Opt. Express 16(6), 3776–3785 (2008). [CrossRef]

19. J. Shah, S. Park, S. Aglyamov, T. Larson, L. Ma, K. Sokolov, K. Johnston, T. Milner, and S. Y. Emelianov, “Photoacoustic imaging and temperature measurement for photothermal cancer therapy,” J. Biomed. Opt. 13(3), 034024 (2008). [CrossRef]

20. P. Saccomandi, E. Schena, and S. Silvestri, “Techniques for temperature monitoring during laser-induced thermotherapy: An overview,” Int. J. Hyperthermia 29(7), 609–619 (2013). [CrossRef]

21. R. Medvid, A. Ruiz, R. J. Komotar, J. R. Jagid, M. E. Ivan, R. M. Quencer, and M. B. Desai, “Current Applications of MRI-Guided Laser Interstitial Thermal Therapy in the Treatment of Brain Neoplasms and Epilepsy: A Radiologic and Neurosurgical Overview,” Am. J. Neuroradiol. 36(11), 1998–2006 (2015). [CrossRef]

22. Z. Tovar-Spinoza, D. Carter, D. Ferrone, Y. Eksioglu, and S. Huckins, “The use of MRI-guided laser-induced thermal ablation for epilepsy,” Childs Nerv. Syst. 29(11), 2089–2094 (2013). [CrossRef]

23. W. Drexler and J. G. Fujimoto, eds., Optical Coherence Tomography: Technology and Applications (Springer, Cham, 2015), 2nd ed. OCLC: 919429452.

24. B. J. Vakoc, M. Shishko, S. H. Yun, W.-Y. Oh, M. J. Suter, A. E. Desjardins, J. A. Evans, N. S. Nishioka, G. J. Tearney, and B. E. Bouma, “Comprehensive esophageal microscopy by using optical frequency-domain imaging (with video),” Gastrointest. Endoscopy 65(6), 898–905 (2007). [CrossRef]

25. M. J. Suter, M. J. Gora, G. Y. Lauwers, T. Arnason, J. Sauk, K. A. Gallagher, L. Kava, K. M. Tan, A. R. Soomro, T. P. Gallagher, J. A. Gardecki, B. E. Bouma, M. Rosenberg, N. S. Nishioka, and G. J. Tearney, “Esophageal-guided biopsy with volumetric laser endomicroscopy and laser cautery marking: a pilot clinical study,” Gastrointest. Endoscopy 79(6), 886–896 (2014). [CrossRef]

26. M. J. Suter, B. J. Vakoc, P. S. Yachimski, M. Shishkov, G. Y. Lauwers, M. Mino-Kenudson, B. E. Bouma, N. S. Nishioka, and G. J. Tearney, “Comprehensive microscopy of the esophagus in human patients with optical frequency domain imaging,” Gastrointest. Endoscopy 68(4), 745–753 (2008). [CrossRef]

27. K. Kurokawa, S. Makita, Y.-J. Hong, and Y. Yasuno, “Two-dimensional micro-displacement measurement for laser coagulation using optical coherence tomography,” Biomed. Opt. Express 6(1), 170–190 (2015). [CrossRef]

28. W. Chang, Y. Fan, X. Zhang, and H. Liao, “An Intelligent Theranostics Method Using Optical Coherence Tomography Guided Automatic Laser Ablation for Neurosurgery,” in 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), (2018), pp. 3224–3227. ISSN: 1558-4615, 1557-170X.

29. W. C. Y. Lo, N. Uribe-Patarroyo, A. S. Nam, M. Villiger, B. J. Vakoc, and B. E. Bouma, “Laser thermal therapy monitoring using complex differential variance in optical coherence tomography,” J. Biophotonics 10(1), 84–91 (2017). [CrossRef]

30. W. C. Y. Lo, N. Uribe-Patarroyo, K. Hoebel, K. Beaudette, M. Villiger, N. S. Nishioka, B. J. Vakoc, and B. E. Bouma, “Balloon catheter-based radiofrequency ablation monitoring in porcine esophagus using optical coherence tomography,” Biomed. Opt. Express 10(4), 2067–2089 (2019). [CrossRef]

31. C. Lee, G. Cheon, D.-H. Kim, and J. U. Kang, “Feasibility study: protein denaturation and coagulation monitoring with speckle variance optical coherence tomography,” J. Biomed. Opt. 21(12), 125004 (2016). [CrossRef]

32. B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Real-time microscopic visualization of tissue response to laser thermal therapy,” J. Biomed. Opt. 12(2), 020501 (2007). [CrossRef]

33. T. J. Pfefer, B. Choi, G. Vargas, K. M. McNally-Heintzelman, and A. J. Welch, “Mechanisms of laser-induced thermal coagulation of whole blood in vitro,” in Lasers in Surgery: Advanced Characterization, Therapeutics, and Systems IX, vol. 3590 (International Society for Optics and Photonics, 1999), pp. 20–31.

34. P. Gong, S. Es’haghian, K.-A. Harms, A. Murray, S. Rea, B. F. Kennedy, F. M. Wood, D. D. Sampson, and R. A. McLaughlin, “Optical coherence tomography for longitudinal monitoring of vasculature in scars treated with laser fractionation,” J. Biophotonics 9(6), 626–636 (2016). [CrossRef]

35. J. F. d. Boer, S. M. Srinivas, A. Malekafzali, Z. Chen, and J. S. Nelson, “Imaging thermally damaged tissue by polarization sensitive optical coherence tomography,” Opt. Express 3(6), 212–218 (1998). [CrossRef]

36. K. Schoenenberger, B. W. Colston, D. J. Maitland, L. B. D. Silva, and M. J. Everett, “Mapping of birefringence and thermal damage in tissue by use of polarization-sensitive optical coherence tomography,” Appl. Opt. 37(25), 6026–6036 (1998). [CrossRef]

37. X. Fu, Z. Wang, H. Wang, Y. T. Wang, M. W. Jenkins, and A. M. Rollins, “Fiber-optic catheter-based polarization-sensitive OCT for radio-frequency ablation monitoring,” Opt. Lett. 39(17), 5066–5069 (2014). [CrossRef]

38. I. Ahmad, A. Gribble, I. Murtza, M. Ikram, M. Pop, and A. Vitkin, “Polarization image segmentation of radiofrequency ablated porcine myocardial tissue,” PLoS One 12(4), e0175173 (2017). [CrossRef]

39. H. H. Müller, L. Ptaszynski, K. Schlott, C. Debbeler, M. Bever, S. Koinzer, R. Birngruber, R. Brinkmann, and G. Hüttmann, “Imaging thermal expansion and retinal tissue changes during photocoagulation by high speed OCT,” Biomed. Opt. Express 3(5), 1025–1046 (2012). [CrossRef]

40. A. S. Nam, I. Chico-Calero, and B. J. Vakoc, “Complex differential variance algorithm for optical coherence tomography angiography,” Biomed. Opt. Express 5(11), 3822 (2014). [CrossRef]

41. K. Beaudette, H. W. Baac, W.-J. Madore, M. Villiger, N. Godbout, B. E. Bouma, and C. Boudoux, “Laser tissue coagulation and concurrent optical coherence tomography through a double-clad fiber coupler,” Biomed. Opt. Express 6(4), 1293 (2015). [CrossRef]

42. K. Beaudette, M. Strupler, J. Ren, B. E. Bouma, and C. Boudoux, “Radiometric model for coaxial single- and multimode optical emission from double-clad fiber,” Appl. Opt. 57(5), 1110 (2018). [CrossRef]

43. S. Lemire-Renaud, M. Rivard, M. Strupler, D. Morneau, F. Verpillat, X. Daxhelet, N. Godbout, and C. Boudoux, “Double-clad fiber coupler for endoscopy,” Opt. Express 18(10), 9755–9764 (2010). [CrossRef]

44. W.-J. Madore, E. D. Montigny, O. Ouellette, S. Lemire-Renaud, M. Leduc, X. Daxhelet, N. Godbout, and C. Boudoux, “Asymmetric double-clad fiber couplers for endoscopy,” Opt. Lett. 38(21), 4514–4517 (2013). [CrossRef]

45. X. Attendu, M. H. Bourget, M. P. d. Sivry-Houle, and C. Boudoux, “Coregistered optical coherence tomography and frequency-encoded multispectral imaging for spectrally sparse color imaging,” J. Biomed. Opt. 25(3), 032008 (2019). [CrossRef]

46. K. Beaudette, N. Godbout, and C. Boudoux, “Advances in Multimodal Imaging Using Double-Clad Fiber Couplers,” J. Lightwave Technol. 37(22), 5674–5685 (2019). [CrossRef]

47. R. Guay-Lord, X. Attendu, K. L. Lurie, L. Majeau, N. Godbout, A. K. Bowden, M. Strupler, and C. Boudoux, “Combined optical coherence tomography and hyperspectral imaging using a double-clad fiber coupler,” J. Biomed. Opt. 21(11), 116008 (2016). [CrossRef]

48. D. Lorenser, B. C. Quirk, M. Auger, W.-J. Madore, R. W. Kirk, N. Godbout, D. D. Sampson, C. Boudoux, and R. A. McLaughlin, “Dual-modality needle probe for combined fluorescence imaging and three-dimensional optical coherence tomography,” Opt. Lett. 38(3), 266–268 (2013). [CrossRef]

49. J. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts & Company, 2007).

50. K. Kurokawa, S. Makita, and Y. Yasuno, “Investigation of Thermal Effects of Photocoagulation on Retinal Tissue Using Fine-Motion-Sensitive Dynamic Optical Coherence Tomography,” PLoS One 11(6), e0156761 (2016). [CrossRef]

51. Y. Bromberg, Y. Lahini, E. Small, and Y. Silberberg, “Hanbury Brown and Twiss interferometry with interacting photons,” Nat. Photonics 4(10), 721–726 (2010). [CrossRef]

52. N. Uribe-Patarroyo, A. L. Post, S. Ruiz-Lopera, D. J. Faber, and B. E. Bouma, “Noise and bias in optical coherence tomography signal decorrelation,” OSA Continuum 3, 709–741 (2020). [CrossRef]

53. N. Uribe-Patarroyo, M. Villiger, and B. E. Bouma, “Quantitative technique for robust and noise-tolerant speed measurements based on speckle decorrelation in optical coherence tomography,” Opt. Express 22(20), 24411 (2014). [CrossRef]

54. S. Makita, K. Kurokawa, Y.-J. Hong, M. Miura, and Y. Yasuno, “Noise-immune complex correlation for optical coherence angiography based on standard and Jones matrix optical coherence tomography,” Biomed. Opt. Express 7(4), 1525–1548 (2016). [CrossRef]

55. R. Shao and L. Chang, “A new maximum power point tracking method for photovoltaic arrays using golden section search algorithm,” in 2008 Canadian Conference on Electrical and Computer Engineering, (2008), pp. 000619–000622. ISSN: 0840–7789.

56. V. Y. Zaitsev, L. A. Matveev, A. L. Matveyev, G. V. Gelikonov, and V. M. Gelikonov, “A model for simulating speckle-pattern evolution based on close to reality procedures used in spectral-domain OCT,” arXiv:1406.3448 [physics] (2014).

57. J. H. Fugett, H. Bennett, J. Shrout, and J. Coad, “Nitro Blue Tetrazolium (NBT) Staining as an Assessment of Thermal Penetration,” The FASEB J. 31, 984 (2017).

58. S. Zhang, “Finding the Right Color Palettes for Data Visualizations,” (2018).

59. L. M. Oliveira, A. L. V. S. Lage, M. P. P. Clemente, and V. V. Tuchin, “Rat muscle opacity decrease due to the osmosis of a simple mixture,” J. Biomed. Opt. 15(5), 055004 (2010). [CrossRef]

60. L. Oliveira, A. Lage, M. Pais Clemente, and V. Tuchin, “Optical characterization and composition of abdominal wall muscle from rat,” Opt. Lasers Eng. 47(6), 667–672 (2009). [CrossRef]

61. P. Kovesi, “Good Colour Maps: How to Design Them,” arXiv:1509.03700 [cs] (2015). ArXiv: 1509.03700.

62. E. Collett, Field Guide to Polarization (SPIE, 2005).

63. N. Uribe-Patarroyo and B. E. Bouma, “Rotational distortion correction in endoscopic optical coherence tomography based on speckle decorrelation,” Opt. Lett. 40(23), 5518 (2015). [CrossRef]

64. A. J. Welch and M. J. V. Gemert, eds., Optical-Thermal Response of Laser-Irradiated Tissue (Springer, 2011), 2nd ed.

65. J. Lee, W. Wu, J. Y. Jiang, B. Zhu, and D. A. Boas, “Dynamic light scattering optical coherence tomography,” Opt. Express 20(20), 22262–22277 (2012). [CrossRef]

66. N. Weiss, T. G. v. Leeuwen, and J. Kalkman, “Localized measurement of longitudinal and transverse flow velocities in colloidal suspensions using optical coherence tomography,” Phys. Rev. E 88(4), 042312 (2013). [CrossRef]

67. N. Uribe-Patarroyo and B. E. Bouma, “Velocity gradients in spatially-resolved laser Doppler flowmetry and dynamic light scattering with confocal and coherence gating,” Phys. Rev. E 94(2), 022604 (2016). [CrossRef]

68. R. Maltais-Tariant, “Real-time OCT surveillance of laser therapy through speckle decorrelation,” Master’s thesis, École Polytechnique de Montréal (2019).

69. D. G. Lowe, “Distinctive Image Features from Scale-Invariant Keypoints,” Int. J. Comput. Vis. 60(2), 91–110 (2004). [CrossRef]

70. B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” Proc. 7th Int. Jt. Conf. on Artif. Intell. (IJCAI) 2, 674–679 (1981).

71. D. Wagner, A. Mulloni, T. Langlotz, and D. Schmalstieg, “Real-time panoramic mapping and tracking on mobile phones,” in 2010 IEEE Virtual Reality Conference (VR), (2010), pp. 211–218.

72. M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33(2), 156 (2008). [CrossRef]

73. H. S. Stone, M. T. Orchard, E.-C. Chang, and S. A. Martucci, “A fast direct Fourier-based algorithm for subpixel registration of images,” IEEE Trans. Geosci. Remote Sensing 39(10), 2235–2243 (2001). [CrossRef]

74. B. Reddy and B. Chatterji, “An FFT-based technique for translation, rotation, and scale-invariant image registration,” IEEE Trans. on Image Process. 5(8), 1266–1271 (1996). [CrossRef]

Name	Description
Visualization 1	Progression of the coagulation highlighted by the correlation of the speckle pattern over 3.5 seconds of exposure by 10~mJ pulses therapy laser at 100~Hz. A label shows when the therapy laser is on. Scale bars are 0.5~mm.
Visualization 2	Highlighted coagulation during laser therapy was automatically stopped when the targeted depth was reached with histology comparison and OCT tomogram for comparison. The red bar shows the targeted depth. A label shows when the therapy laser is on. Sc

	Fast Mode	Retrospective Diagnostic Mode	Depth Targeting Mode
B-scan Frequency	1 kHz	100 Hz	100 Hz
B-scan Width	0.35 mm	3.71 mm	3.71 mm
B-scan shape	Triangular wave	Sawtooth wave	Sawtooth wave
Max Laser Power	CW: 1.66 W	10 mJ at 100 Hz 3 ms pulse width	10 mJ at 100 Hz 3 ms pulse width
Treated Area	Entire lateral field of view	30% of the lateral field of view	30% of the lateral field of view
Laser Shutter	Triggered by timer	Triggered by timer	Triggered when targeted coagulation depth reached
Saved Data	OCT fringe data	OCT fringe data	OCT intensity
Typical $Δ$ z	10 px or 77.9 µm	3 px or 32.5 µm	3 px or 32.5 µm
Typical $Δ$ x	16 px or 237 µm	10 px or 46.9 µm	50 px or 234.4 µm
Typical $Δ$ t	100 px or 100 ms	10 px or 100 ms	10 px or 100 ms

Variable	Value	Description
$p_{x}$	10 µm	pixel physical size
$s_{z}$	120	Number of pixels in $z$ direction
$s_{x}$	120	Number of pixels in $x$ direction
$s_{n}$	10	Number of scatterers per pixel
$m_{a}$	0.1 µm	Maximum amplitude of scatterers random movements
$τ$	0,1,2, $\dots$ ,200 frames	Delay of $g^{(2)}$
$Δ$	0,0.025,0.05, $\dots$ ,5 pixels	Global movement amplitude
$Δ_{max}$	5 pixels	Maximum global movement amplitude
$w$	4.25 µm	PSF beam waist

	Fast Mode	Retrospective Diagnostic Mode	Depth Targeting Mode
B-scan Frequency	1 kHz	100 Hz	100 Hz
B-scan Width	0.35 mm	3.71 mm	3.71 mm
B-scan shape	Triangular wave	Sawtooth wave	Sawtooth wave
Max Laser Power	CW: 1.66 W	10 mJ at 100 Hz 3 ms pulse width	10 mJ at 100 Hz 3 ms pulse width
Treated Area	Entire lateral field of view	30% of the lateral field of view	30% of the lateral field of view
Laser Shutter	Triggered by timer	Triggered by timer	Triggered when targeted coagulation depth reached
Saved Data	OCT fringe data	OCT fringe data	OCT intensity
Typical $Δ$ z	10 px or 77.9 µm	3 px or 32.5 µm	3 px or 32.5 µm
Typical $Δ$ x	16 px or 237 µm	10 px or 46.9 µm	50 px or 234.4 µm
Typical $Δ$ t	100 px or 100 ms	10 px or 100 ms	10 px or 100 ms

Variable	Value	Description
$p_{x}$	10 µm	pixel physical size
$s_{z}$	120	Number of pixels in $z$ direction
$s_{x}$	120	Number of pixels in $x$ direction
$s_{n}$	10	Number of scatterers per pixel
$m_{a}$	0.1 µm	Maximum amplitude of scatterers random movements
$τ$	0,1,2, $\dots$ ,200 frames	Delay of $g^{(2)}$
$Δ$	0,0.025,0.05, $\dots$ ,5 pixels	Global movement amplitude
$Δ_{max}$	5 pixels	Maximum global movement amplitude
$w$	4.25 µm	PSF beam waist

Real-time co-localized OCT surveillance of laser therapy using motion corrected speckle decorrelation

Abstract

1. Introduction

2. Theoretical model

2.1 Speckle decorrelation during coagulation

2.2 Noise correction

2.3 Motion correction

3. Materials and methods

3.1 Experimental setup

3.2 Software

3.3 Histology

4. Results

4.1 Noise and motion correction

4.2 Distinguishing coagulated from coagulating tissues

4.3 Monitoring in real-time

5. Discussion

6. Conclusion

7. Appendix A - Speckle simulation

Funding

Acknowledgments

Disclosures

References

Supplementary Material (2)

Cited By

Figures (15)

Tables (2)

Equations (23)

Biomedical Optics Express