Analysis and evaluation of BC-mode OCT image visualization for microsurgery guidance

: Optical coherence tomography (OCT) has been gaining acceptance in image-guided microsurgery as a noninvasive imaging technique. However, when using B-mode OCT imaging, it is diﬃcult to continuously keep the surgical tool in the imaging ﬁeld, and the image of the tissue beneath the tool is corrupted by shadow eﬀects. The alternative using C-mode OCT imaging is either too slow in imaging speed when operating in a high-resolution mode, or provides a poor image resolution in a high-speed mode, with the sweep rate less than one million hertz. Moreover, the 3-dimensional rendering of C-mode OCT image makes it diﬃcult to visualize the tissue structure and track the surgical tool beneath the tissue surface. To solve these problems, we propose a BC-mode OCT image visualization method. This method uses a sparse C-scanning scheme, which provides a set of high-resolution B-mode OCT images at sparsely spaced cross sections. The ﬁnal BC-mode OCT image is obtained by averaging the image set, with inter frame variance processing to enhance the signal of the surgical tool and tissue layers. The performance of BC-mode OCT images, such as image resolution, signal to noise ratio (SNR), imaging speed, and surgical tool tracking accuracy, is analyzed theoretically and veriﬁed experimentally. The feasibility of the proposed method is evaluated by guiding the insertion of a 30-gauge needle into the cornea of an ex-vivo human eye freehand. The results show that this provides better visualization of both the surgical tool and the tissue structure than the conventional B- or C-mode OCT image.


Introduction
Microsurgery, such as ophthalmic surgery and neurosurgery, involves manipulation of delicate tissues and requires precise control of the surgical tool [1]. Therefore, providing a clear visualization of both the surgical tool and the tissue structure is of great importance to the microsurgery guidance. Conventionally, surgeons use surgical microscopes to magnify the field of view, but the view is restricted to an en face perspective with limited depth information [2]. Alternatively, some tomographic imaging modalities, such as ultrasound, computed tomography (CT) and magnetic resonance imaging (MRI) provide structure information inside the tissue but still suffer from limited resolution [3], radiation exposure [4] and incompatibility with surgical tools [5].
Optical coherence tomography (OCT) [6], a high-speed imaging modality with micrometer resolution, has been gaining acceptance in image-guided microsurgery [7][8][9]. Different visualization modes for OCT images have been explored to better guide the surgical tool. Kang et al. used the A-mode common path OCT image to guide the microinjector [10]. Shin et al. used the M-mode OCT image to guide the surgical needle for lamellar keratoplasty [11]. Yu et al. used B-mode OCT image for robotic microsurgery [12]. Kang et al. [7], Viehland et al. [13] and Keller et al. [14] used a C-mode OCT image for surgical tool tracking. This work shows that A-mode and M-mode OCT images only provide depth information, and thus their visualization can be difficult to interpret. C-mode OCT images are easier to interpret. However, with the A-scan rate less than one million hertz, C-mode OCT imaging is either slow when the imaging speed is less than 5 Hz, or poor in image resolution when the imaging speed is higher than 5 Hz, which is necessary in order to provide a decent field of view [15][16][17][18][19]. Although video rate C-mode OCT imaging [20] and densely sampled wide-field snapshot C-mode OCT imaging [21] have been demonstrated, it is at the cost of using an ultra high-speed swept source on the order of million hertz. It is not fair to compare when the sweep rate is different. Therefore, a better criteria is the scanning spots within each frame. Usually, compared with B-mode OCT images, C-mode OCT images need far more scanning spots within each frame. Therefore, B-mode OCT imaging is usually much faster than C-mode OCT image, under the condition of the same sweep rate. Moreover, without projection to a 2-dimensional plane, the 3-dimensional rendering of C-mode OCT image highlights the tissue surface, so it is difficult to visualize the tissue structure and track the surgical tool beneath the surface. Compared with the C-mode OCT image, B-mode OCT image provides a clearer cross-sectional view of tissue structure and surgical tool underneath the tissue surface, and it is faster in imaging speed. Therefore, the B-mode OCT image is a better intraoperative imaging mode than C-mode OCT image. However, in conventional B-mode OCT imaging, it is difficult to continuously keep the surgical tool in the imaging field since it provides only one cross-sectional image of a single plane. Moreover, when the tool is in the imaging plane, the tissue signal is blocked by the surgical tool, which results in shadow effects in the B-mode OCT image.
In this paper, we propose a BC-mode OCT image visualization method to solve these problems, which compresses the sparsely scanned OCT volume data into a compounded B-mode OCT image. First, we introduce the BC-mode OCT system setup, scanning configuration, and calculation method. Then, we analyze the performance of the BC-mode OCT imaging theoretically and verify the theoretical results with experimental data. Next, we describe in detail more image processing techniques that we use to enhance the BC-mode OCT image quality. Finally, we evaluate the BC-mode OCT image performance by guiding a 30-gauge needle into the ex-vivo human cornea freehand. The results show that this provides better visualization than the conventional Bor C-mode OCT image, which validates its effectiveness in guiding the surgical tool during microsurgery.

System setup
The system setup is shown in Fig. 1(a). The system uses a 100kHz swept source OCT centered at 1060nm with 110nm tuning range. The output of the swept laser is split into two arms by a 75:25 fiber coupler. The light source couples to the sample arm with 25%, and to the reference arm with 75%. The output optical power in the sample arm is 2.5mW. Polarization controllers are used to maintain the same light polarization states between the two arms. Achromatic collimators are used to collimate the light from the fiber. Galvanometer scanners, which contain two scanning mirrors (shown in Fig. 1(b)), are used for scanning. The objective lens is used to focus the collimated light onto the sample. The numerical aperture (NA) of the objective lens is 0.05. The dispersion compensation lens is used to compensate for the dispersion due to the light double passing through the object lens. The signal from both the sample arm and the reference arm are coupled by a 50:50 fiber coupler and detected by a balanced detector. The interferometric data is processed by a data acquisition board (12-bit, 500 MSPS) and then transferred to the workstation by a frame grabber for further data processing. Graphics processing units (GPU) are used for high-speed data processing [16,22]. The system is controlled by our customized software programmed by C++ and C#.
We use a standard OCT C-scanning, which is shown in Fig. 1(b). Two scanning mirrors are used for scanning. The first one scans densely along the X axis, while the second one scans sparsely along the Z axis. The orientation and the movement of the surgical tool are both set to be roughly parallel to the densely scanning axis X. The yellow lines shown in Fig. 1 the multiple scanning cross sections in the XY plane. After one scanning period, the system can obtain a set of OCT volume data, which is an image set that contains multiple B-mode OCT images from different cross sections.
The BC-mode OCT image is obtained by averaging these multiple B-mode OCT images from different cross sections, which is calculated by Eq. (1), where I BC is the BC-mode OCT image, S k is the k-th image in the image set, N is the total number of images in the image set and (i, j) indicates the pixel index in the image. Although averaging techniques have been widely used in OCT imaging, they are usually performed by taking the average of multiple scans at the same cross section [23,24], while BC-mode averaging scheme is taking the average from different cross sections. Some similar averaged B-scan images from different cross sections were previously demonstrated by Tao et al. [25] and Ehlers et al. [26], but without theoretical analysis or further discussion about its performance. Moreover, the results still suffer from the tool shadow effects due to the narrow scanning range. To solve this problem, we made the scanning along the Z axis sparse enough to cover a larger scanning range. The specific scanning configuration and the rationale are discussed in detail in the following paragraph.
The scanning configuration is shown in Table 1. The illustration of the scanning region and each labeled axis XYZ is shown in Fig. 1(b). Along the X axis, the scanning length is 6000 µm with a 6 µm dense scanning step. Along the Z axis, the scanning length is 2000 µm with a 250 µm sparse scanning step. Therefore, a set of OCT volume data contains 8 B-mode OCT images in total. The B-scan rate is 100Hz and the C-scan rate is 12.5Hz. Since one BC-mode OCT image is calculated from one set of OCT volume data, the final BC-mode OCT imaging speed is 12.5Hz. The size of the OCT volume data is 1000*1024*8 (X*Y*Z). This scanning configuration is specifically customized for the deep anterior lamellar keratoplasty (DALK) procedure, where we need to use our OCT system to guide a 30-gauge needle to Descemet's membrane of the cornea of an ex-vivo human eye. This is also where we collected the data in this paper. To further explain the scanning configuration of our method, we show two adjacent B-mode OCT images from the image set which contains the needle signal in Fig. 2(a) and 2(b). As we can see in both images, the needle blocks out the signal from the tissue beneath it, creating a shadow in the image. Moreover, due to the limit of a single imaging plane, parts of the needle are missing, but then appear in the adjacent B-mode OCT image. Although separately both images show the tissue shadowing effect and the difficulties of keeping the needle in the imaging field, together they provide more complete information about the needle and the cornea structure. The scanning step along Y axis is set to be a little less than the outer diameter of the 30-gauge needle (311 µm) so that the whole needle structure can be captured by the OCT volume data as shown in Fig. 2(a) and 2(b). Moreover, the 6000 µm by 2000 µm total scanning range provides a wide field of view, making it easier to keep the needle position inside the scanning range during the free hand operation. The main idea of BC-mode OCT image visualization is to compress the sparsely scanned OCT volume data into one compounded B-mode OCT image. To understand its advantage and limitations, it is important to compare its performance against the conventional B-mode OCT image visualization. In the following section, we will derive a theoretical model for BC-mode OCT image visualization. Based on this model, we will analyze its performance in terms of the resolution, signal to noise ratio (SNR), imaging speed and needle tracking accuracy.

Resolution
OCT image resolution includes both lateral resolution and axial resolution. For the conventional B-mode OCT image visualization, both axial and lateral resolutions are independent and depend totally on the system parameters. Specifically, when the transversal sampling rate is sufficient, the lateral resolution depends on the diffraction limit of the OCT system [27]. At the focal plane, the lateral resolution can be calculated by Eq. (2), where σ x is the lateral resolution, λ 0 is the center wavelength of the laser, NA is the numerical aperture of the objective lens (shown in Fig. 1(a)). Although the lateral resolution depends on the imaging depth, it does not affect the following analysis on image resolution since its value is not specified in the analysis. The axial resolution depends on the coherence length of the laser. A common definition of the axial resolution is the full-width half maximum (FWHM) of the point spread function (PSF) in the time domain [27,28]. The axial resolution can be calculated by Eq. (3) if we consider the gaussian profile laser spectrum, where σ y is the axial resolution, λ 0 is the center wavelength of the laser, ∆λ is the FWHM of the laser spectrum in λ space. In our OCT system, however, the swept source produces a nearly rectangular spectrum. The constant factor at the right hand side of Eq. (3) will be slightly different. Fourier transformation is used to obtain the time domain PSF. The coherence length is defined as the FWHM of the PSF, while the axial resolution is given by half the coherence length since the wave travels twice the sample arm [29]. Therefore, the constant factor should be modified by 0.6 for rectangular source spectrum. After we have both σ x and σ y , the point spread function for B-mode OCT image system (B-PSF) can be described by Eq. (4), where h B is B-PSF, σ x and σ y are lateral and axial resolution respectively. For the BC-mode images, on the other hand, the image resolution not only depends on the system parameters but also depends on the sample characteristics. The idea of BC-mode OCT image is to average laterally shifted B-mode OCT images. However, it is not ensured that each laterally shifted B-mode OCT image has the exact same structure, in which case the resolution will remain the same. Qualitatively, if the sample has highly flat structures, the influence of the sample structure on the image resolution can be neglected. On the contrary, if the sample structure has a large curvature, it will affect the image resolution significantly.
To quantitatively analyze how the sample structure affects the image resolution, the strategy we use here is to convert the discrete problem into a continuous function space problem. We first define a continuous OCT image function space Ω = {f |f : R 3 → R}, where its element f maps the 3-dimensional continuous position space to the 1-dimensional continuous image intensity space. Then, BC-mode OCT images can be calculated by Eq. (5), where I BC is the BC-mode OCT image, t is the sparsely scanning length along the Z axis, z 0 is the center position of the scanning range along the Z axis, and (x, y, z 0 + z) indicates the position in the 3-dimensional continuous position space. This equation is the continuous version of Eq. (1). Specially, the B-mode OCT image at the z = z 0 cross section can be calculated by taking the limit t → 0 if we neglect the laser beam width. To proceed, we approximate the function within a small integration range by using the first order Taylor expansion, which is shown in Eq. (6), The next step is to find the point spread function for the BC-mode OCT image (BC-PSF), which is a function of x and y. In order to do that, we need to make Eq. (6) integrate on the XY plane instead of integrating along the Z axis. This can be done by finding the contour line with the same image intensity, i.e. x(z), y(z) as a function of z so that Eq. (7) is satisfied, By taking the derivative with respect to z on both side of Eq. (7) and using the chain rule, we can get a differential geometric relationship of the contour line shown in Eq. (8), By plugging Eq. (8) into Eq. (6) and use the inverse Taylor expansion, we can obtain Eq. (9), To analyze the contour line quantitatively, we define two angles α and β, which are shown in Fig. 3. Note that we rotate the axis in order to better illustrate the angle relationship. ì n is the unit tangential vector of the contour line. Angle β is the polar angle between ì n and the Z axis. Angle α is the azimuth angle between the orthogonal projection of ì n on the XY plane and the X axis. Therefore, the unit tangential vector ì n can be described by Eq. (10), ì n = (sin β cos α, sin β sin α, cos β). Then, we can get the differential geometric relationship in Eq. (11) and Eq. (12), Plugging Eq. (11) and Eq. (12) into Eq. (9), the BC-mode OCT image term becomes Eq. (13), Substituting f in Eq. (13) with B-PSF h B in Eq. (4), we can obtain BC-PSF in Eq. (14), where we do not write out the third argument z since it is a constant. The geometric meaning of the polar angle β is the curvature angle of the sample structure along the sparsely scanning axis Z. A smaller β will lead to a smaller effect on the resolution of the BC-mode OCT image. As we can see from Eq. (14), if β = 0, BC-PSF is the same as B-PSF and the image resolution will not be affected. The angle α denotes the direction of the orthogonal projection of the tangential vector on the XY plane. Two extreme situations happen when α = 0 and α = π/2. As we can see from Eq. (14), if α = 0, then sin α = 0 and the axial resolution will not be affected. If α = π/2, then cos α = 0 and the lateral resolution will not be affected. These results are reasonable and expected. For a fixed point in the OCT volume data, the sample structure is fixed, and so are the angles α and β. For fixed angles α and β, Eq. (14) is a line integration. The integration can be rewritten as Eq. (15), where Rect(2z/t) is a rectangular function, i.e. Rect(2z/t) = 1 for −t/2<z<t/2, otherwise Rect(2z/t) = 0. To make Eq. (15) integrable, we approximate the rectangular function as a Gaussian function with the same variance. The reason we do this approximation is that FWHM of a PSF is more important than the shape of the PSF, which determines the image resolution directly. Therefore, instead of numerically obtaining the exact shape of the BC-PSF, we are trying to find the broadened width when the rectangular function convolutes with the PSF based on a statistical approach. By treating the rectangular function as a probability distribution p(x) for a random variable X, we can get the variance of X by V(X) = E(X 2 ) − E(X) 2 = t 2 /12. Then, we approximate the rectangular function in Eq. (16), where σ 2 t = t 2 /12. The integration result for BC-PSF in Eq. (15) is approximated in Eq. (17), where σ 2 0 = σ 2 x σ 2 y /(8 ln 2 · tan 2 β(σ 2 x sin 2 α + σ 2 y cos 2 α)), and the matrix A is shown in Eq. (18), Therefore, finding the image resolution is equivalent to finding the eigenvalue of the matrix A. Note that n x = (1, 0) and n y = (0, 1) are no longer the eigenvectors of the matrix A if β 0 and α 0 or π/2. This is expected since the line integration is not along the axis X or axis Y in this case. However, the eigenvectors are not simply parallel or perpendicular to the integration direction due to the asymmetry along the X axis and the Y axis, i.e. σ x σ y . Despite its complexity, the eigenvalues of the matrix A can still be found but it is not necessary to show the results here. Instead, we can do approximation because for a large scanning length t, the difference between σ x and σ y can be neglected. Thus, we approximate that σ 2 In this approximation, we get symmetry along the X axis and the Y axis, thus the eigen problem becomes much easier to solve. The eigenvalues and the eigenvectors are shown in Eq. (19) and Eq. (20), where λ 1 and λ 2 are the eigenvalues, ì n 1 and ì n 2 are the corresponding eigenvectors, which are parallel or perpendicular to the integration direction. Therefore, along the direction ì n 1 = (cos α, sin α), the resolution term can be found as in Eq. (21), Since ì n 1 is parallel to the integration direction, the resolution should be the worst, and thus it is denoted by σ max . Along the direction ì n 2 = (− sin α, cos α), the resolution term can be found as in Eq. (22), Since ì n 2 is perpendicular to the integration direction, the resolution should be the best, and thus it is denoted by σ min . For any other direction, the image resolution lies between σ min and σ max , i.e. σ min <σ<σ max . This result is important since it approximates the upper bound and lower bound for the image resolution along any direction.
However, there are a lot of points inside the OCT image volume data; each point has its unique sample structure, i.e. the azimuth angle α and the polar angle β. It is also important to analyze the average effect of the sample structure on the image resolution. To do that, we made a statistical assumption that for a fixed polar angle β, the azimuth angle α is distributed uniformly between 0 to π. Therefore, for a fixed angle β, we can average the BC-PSF in Eq. (14) for different angle α as in Eq. (23), where we use the area measure zdαdz and 4/πt 2 is the normalization factor. We change the integrate range for α from [0, π] to [0, 2π] because we change the integration range for z from After transforming the basis from the polar coordinate to the Cartesian coordinate, where x 0 y 0 = z cos α z sin α , the integration can be written as Eq. (24), where Circle(( To make the Eq. (24) integrable, we again approximate the circular function as a 2-dimensional Gaussian function. This time, we regard the circular function as a joint probability distribution p(x, y) for random variables X and Y. To find the variance for random variable X, we need to calculate the marginal probability distribution by p X (x) = ∫ +∞ −∞ dyp(x, y). Intuitively, the 2-dimensional circular function makes X more confined to the center than the 1-dimensional rectangular function, thus it should have smeller variance. We can calculate the variance of X by V(X) = E(X 2 ) − E(X) 2 = t 2 /64, which is smaller than t 2 /12 as expected. Therefore, we approximate the circular function in Eq. (25), 4 πt 2 Circle where σ 2 t = t 2 /64. The integration result for BC-PSF in Eq. (24) is approximated in Eq. (26), Therefore, the average effects of the sample structure on the lateral image resolution and axial image resolution are separated and can be described by Eq. (27) and Eq. (28), whereσ x andσ y lie within σ min and σ max for large t as expected. Note that although the lateral resolution σ x depends on the imaging depth, Eq. (21), Eq. (22), Eq. (27) and Eq. (28) still show the correct trend for the BC-mode OCT image resolution compared with the conventional B-mode OCT image at the same imaging depth.
To validate our theory, we used the surface of the cornea as landmarks to approximate the BC-mode OCT image resolution. The human cornea surface usually has a spherical shape with a 6.75 mm radius. The sparsely scanning length along the axis Z is 2 mm in our scanning configuration. Figure 4 shows the model of the human cornea, where r is the radius of the cornea and t is the scanning length along the Z axis. Considering that we are imaging the center part of the cornea, the center of the cornea is aligned with the center of the scanning region along the Z axis. The minimal polar angle of the cornea structure is zero at the center of the scanning region, while the maximal polar angle is tan −1 (t/2r) at the edge of the scanning region. Statistically, we want to obtain an average polar angle, which is denoted by the blue dashed line and angle β in Fig. 4. Since this polar angle is small in cornea structure, we can approximate the average polar angle β by Eq. (29), Therefore, we get the approximate average polar angle β = 4.2 • . For the OCT system parameters, we use the center wavelength λ 0 = 1060nm, the wavelength tuning range ∆λ = 110nm, the numerical aperture of the scanning lens NA = 0.05. We can approximate the lateral resolution by using Eq. (2) since the imaging plane is close to the focal plane of the laser beam. We tested three different regions of interest (ROIs) around the center part of the cornea, which are shown in Fig. 5(a)-5(c). Figure 5  To calculate the resolution of the BC-mode OCT image, we segmented the surface of the cornea in the B-mode OCT image. The segmentation method we use here is based on graph theory and dynamic programming, with the accuracy around ±1 pixel [14,30]. The resolution for each A-line is approximated by calculating the largest difference between the segmented surface positions. To get the dependence of the resolution on the sparsely scanning length, we acquire BC-mode OCT images for different N from 2 to 8 so that the sparse scanning length ranges from 0.5 mm to 2 mm. In Fig. 5(d)-5(f), the black circles show the mean value of the resolution, while the black bars show the variance value of the resolution. We can see that most resolution values lie within the lower bound and upper bound. Moreover, for each tested data, the resolution has a nearly linear relationship with the sparsely scanning length. The difference in the slope comes from another structure parameter, the azimuth angle α. As we discussed before, when α is close to 0, the resolution is close to the lower bound; when α is close to π/2, the resolution is close to the upper bound. Therefore, we can say that ROI-1 has a lower value of α, while ROI-2 has a larger value of α.

SNR
One benefit of BC-mode averaging is to reduce the speckle noise. Since the signal to noise ratio (SNR) is a good quantification for the level of speckle noise, we need to compare the SNR between BC-mode OCT images and conventional B-mode OCT images. First, we define the SNR as shown in Eq. (30), where µ s is the mean intensity of the signal, σ n is the standard deviation of the noise. Since the BC-mode OCT image is averaged from many laterally shifted B-mode OCT images, the mean intensity µ BC of the signal in the BC-mode image should remain the same as the mean intensity To prove this, we can use the continuous OCT image function space Ω that we defined in the resolution section. The mean intensity µ B of B-mode OCT image can be calculated in Eq. (31), where we integrate over the signal range and S is the total integration area, f ∈ Ω, and the plane z = z 0 is the scanning cross section of the B-mode OCT image. By using the continuous BC-mode OCT integration scheme in Eq. (5), the mean intensity µ BC of the signal in the B-mode image can be calculated in Eq. (32), where z = z 0 is the center plane of the scanning range along the Z axis, and t is the sparsely scanning length along the Z axis. This is true for BC-mode OCT signals that are not blocked by the surgical tool. However, for BC-mode OCT signals that are blocked by the surgical tool, we need to modify Eq. (32), which we will discuss explicitly in the later part of this section. By exchanging the integration order and plugging Eq. (31), we can get Eq. (33), To proceed, an assumption was made that for a small range of area, the mean intensity of the B-mode signal remains the same, i.e. µ B (z 0 + z) = µ B (z 0 ) for any small z. Under this assumption, we can easily get Eq. (34), where we prove that the mean intensity µ BC should remain the same as µ B and it does not depend on the averaging frame number N. However, the case for the standard deviation of the noise σ n is different and it will not remain the same for different N. In fact, OCT speckle noise is spatially independent, and this independence is widely used in speckle reduction of OCT image [31][32][33]. Therefore, we assume that the noises for different B-mode OCT images X 1 , X 2 · · · X N are random variables from an identical independent distribution with a standard deviation σ B , i.e. X 1 , X 2 · · · X N ∈ i.i.d and V(X i ) = σ B for i = 1, 2, · · · , N, where N is the number of averaging frames. Then the noise in a BC-mode OCT image can be modelled as Eq. (35),X The variance of the noiseX in BC-mode OCT image can be calculated in Eq. (36), We label the image area where OCT signals are not blocked by the surgical tool the unshaded area. By using Eq. (30), Eq. (34) and Eq. (36), the SNR of the unshaded area can be modeled as Eq. (37), where SNR B = 10log 10 (µ 2 B /σ 2 B ) is the SNR of the conventional B-mode OCT image. However, for the image area where OCT signals are blocked by the surgical tool, which we label the shaded area, the SNR is different and should be lower than that of the unshaded area. Since the shaded area is filled with noise instead of the OCT signal, Eq. (33) and Eq. (34) do not work. To model it, we first define t 0 as the thickness of the surgical tool along the sparsely scanning axis Z. Therefore, Eq. (33) can be easily modified as Eq. (38), where µ n is the mean intensity of the noise at the shaded area, and −(t − t 0 )/2<a<(t − t 0 )/2 since the surgical tool should be kept within the sparsely scanning range. We can ignore the independent variables z 0 + z based on the assumption that both µ B and µ n remain constant for a small range of area. Therefore, the mean intensity µ BC of signal in the BC-mode OCT image can be written in Eq. (39), We define another parameter s as the scanning step so that t = Ns. Therefore, the SNR of the shaded area can be modeled as Eq. (40), It is generally true that in OCT image, the signal level is higher than the noise level, i.e. µ B >µ n . Therefore, (µ B − µ n )/µ B >0 and we can see that the SNR increases when the scanning step s increases. However, in practice, the scanning size s should be set less than the thickness of the surgical tool t 0 , i.e. s ≤ t 0 . Therefore, when s = t 0 , the best SNR is obtained for a shaded area which is shown in Eq. (41), where we can see that the noise in the shaded area µ n is not a bad thing since it improves the SNR. If the mean intensity of the noise is the same as the mean intensity of the signal, i.e. µ n = µ B , then Eq. (11) will go back to Eq. (9) as expected. Usually the mean intensity of the noise µ n changes with the axial position and it is smaller at a deeper axial position. However, no matter how µ n changes, it should be bounded by 0 and µ B , i.e. 0<µ n <µ B . Therefore, the SNR of the shaded area is also bounded. The upper bound of the SNR in the shaded area is given by Eq. (42) when µ n = µ B , SNR BC−shaded−max = SNR B + 10log 10 N, which is the same as the SNR of the unshaded area. The lower bound of the SNR is given by Eq. (43) when µ n = 0, which we are more interested in because this defines the worst case of how the shading effect of the surgical tool impacts on the SNR of BC-mode OCT images. To validate our theory on SNR, we compared both theoretical SNR predictions and experimental SNR values calculated from BC-mode images. To obtain the dependence of SNR on the averaging frame number N, we acquire BC-mode OCT images for different N from 1 to 8. As shown in Fig. 6(a), we choose three different ROIs on the BC-mode image. The ROI-1 is used to calculate the standard deviation of the noise. The ROI-2 is used to calculate the mean intensity of the unshaded area, while the ROI-3 is used to calculate the mean intensity of the shaded area. For better comparison, we chose the ROI-2 and the ROI-3 at the same axial depth, as we can see in Fig. 6(a). Figure 6(a) was rescaled to better arrange the figure. To explore the lower bound of the SNR for the shaded area, we made the noise intensity of the shaded area zero. The results are plotted in Fig. 6(b). The blue curve and the blue circles show the theoretical and experimental value for SNR of the unshaded area, respectively. The red curve and the red circles show the theoretical lower bound and experimental value for SNR of the shaded area, respectively. The constant SNR B was found to be 9.43 by least square fitting our theoretical models with the experimental data. Note that the absolute value of SNR is not as important as the relative trend of SNR, since the absolute value of SNR can be easily adjusted by changing the reference arm power or by the background subtraction, while the trend of SNR remains the same. As we can see, the theoretical value and the experimental value match well. The little offset between the theoretical and experimental values is because of the mean intensity µ B of the OCT signal does not remain exactly the same in a small range around the ROIs. From the plot in Fig. 6(b), we can see that when we increase the averaging frame number N for BC-mode OCT image, the SNR for both unshaded and shaded area increases and the difference between the two SNR decreases. However, there are other aspects that we should consider when choosing a proper value of N, such as the imaging resolution as we discussed before, and the imaging speed that we will discuss in the following section.

Imaging speed
The imaging speed for the BC-mode OCT image visualization is easy to analyze. However, it is necessary to be analyzed since it gives us an upper bound for the number of averaging frames N if we want a decent video frame rate. The frame rate is determined by the sweep rate of the laser source and the number of the scanning spots within one frame, which can be described by Eq. (44), where R A is the available A-scan sweep rate of the swept source, and N BC is the scanning spots within one BC-mode OCT image frame. In our OCT system, the swept source has a 100kHz sweep rate, i.e. R A = 100kHz. Since each B-mode OCT image contains 1000 A-scan lines and each BC-mode OCT image contains N B-mode OCT images, we have N BC = 1000 · N. Therefore, by using Eq. (44), the frame rate in our system is 100/N. The frame rate should be at least 10 Hz for guiding the surgical tool at a video rate. Therefore, the averaging number of frames for each BC-mode OCT image should be no more than 10, i.e. N ≤ 10.

Tracking accuracy
Although BC-mode visualization can improve the visual experience by continuously keeping the needle position in the imaging field, it is unknown at what accuracy the needle can be tracked by using BC-mode OCT images. For procedures involving the needle guiding, the most important target requiring precise tracking is the needle tip so that the distance between the needle tip and the target surface such as Descemet's membrane in the case of DALK can be accurately calculated. To quantitatively analyze the needle tip tracking accuracy, we need to consider the needle shape and the BC-mode scanning scheme. Figure 7 shows the geometric model for needle tracking accuracy analysis. In Fig. 7(a), the blue part represents the needle and the red lines indicate the sparse BC-mode scanning cross-sections. The needle tip is modeled as isosceles triangle shape for first order approximation. The distance between the adjacent scanning line is the scanning step size s, which was defined previously. The scanning step size s should satisfy that 0<s<t 0 , where t 0 is the thickness of the surgical tool along the sparsely scanning axis Z. The angle θ is the angle of the needle tip seen from the top view. Figure 7(b) shows the 3-dimensional diagram of the needle, where the top needle is the real needle and the bottom needle is the orthogonal projection of the real needle onto the horizontal plane XZ, which is also shown in the top view in Fig. 7(a). The angle α is the real angle of the needle tip and it is a fixed value for a fixed type of needle. The angle β is the insertion angle between the needle and the horizontal plane. The angle θ is the angle of the needle tip after its projection onto the horizontal plane, which is the same as the one shown in Fig. 7(a). Based on the geometric relationship in Fig. 7(b), the angle θ can be calculated as in Eq. (45), During the surgical procedure, the needle can be randomly placed, while the positions of BC-mode scanning lines are fixed. In practice, we try to make both the needle orientation and the needle movement parallel to the BC-mode scanning axis X during the experiment, so we neglect the small angle that might occur between them in the theoretical analysis. We define x as the distance between the needle tip and the closest BC-mode scanning line. Therefore, x should satisfy that 0 ≤ x ≤ s/2. Moreover, we assume that x is uniformly distributed in this range. The next step is to calculate the measurement error of the needle tip position by using the BC-mode OCT image. The vertical blue dashed line at the left side in Fig. 7(a) shows the real position of the needle tip, while the one at the right side shows the measured position of the needle tip in BC-mode OCT image. The error e is the distance between the two blue dashed lines. Since we care more about the vertical distance between the needle tip and the target surface, we are more interested in the vertical component of the error, e v = e tan β. Using the geometric angle relation in Eq. (45), we can get the vertical error e v in Eq. (46), Since the distance from the needle tip to the closest BC-mode scanning line is bounded, i.e. 0 ≤ x ≤ s/2, the vertical error e v is also bounded. This is the advantage of the BC-mode OCT scanning scheme, where it is much easier to keep the needle in the range of the scanning field compared with B-mode OCT scanning scheme. As long as the needle is kept inside the scanning field, the measurement error of the needle tip is always bounded. However, for B-mode OCT scanning scheme, since there is only one scanning line, the measurement error of the needle tip is not bounded and can be comparably large if the needle does not exactly lie in the B-mode scanning plane.
To proceed, based on our assumption that x is a random variable that is uniformly distributed in the range [0, s/2], we can calculate the mean of the vertical error in Eq. (47), We can also calculate the standard deviation of the vertical error in Eq. (48), The results show that both the mean and the variance of the vertical error increase linearly with the scanning step size s, and they also increase when we increase the insertion angle β. Therefore, with a smaller insertion angle, we will get better tracking accuracy. Furthermore, we can see that the error is always positive, which means that the real needle tip position is always ahead of the measured needle position. Thus, we should compensate for that error by adding a small vertical distance d to the measured needle tip position in the BC-mode OCT image. To find the appropriate value for d, we used minimum variance estimation (MVE). The variance of error after adding d can be described by Eq. (49), Therefore, to minimize the variance σ 2 d , we should choose d = µ e v . In this case, the minimum variance of the measurement error will be the same as the variance of the vertical error, i.e. σ 2 d = σ 2 e v . To better understand how the scanning step size s and the insertion angle β affect the tracking accuracy, the measurement error variances after the distance compensation versus the scanning step size s for different insertion angles β are plotted in Fig. 8. The angle of the 30-gauge needle tip α was measured to be about 60 • . As we can see in Fig. 8, the smaller the insertion angle, the better the accuracy for the needle tracking. Moreover, the tracking error of the needle tip in the BC-mode OCT image increases linearly when we increase the scanning step size s. Although we discussed before that larger scanning step size s will give better SNR, we also need to consider the tracking accuracy when choosing a proper value of s.
In conclusion, the BC-mode OCT image has better tracking accuracy than B-mode OCT image. To further improve the tracking accuracy of BC-mode OCT image, one thing we have to do is to compensate for the error by adding µ e v to the actual measured needle tip vertical position, so that we will have minimal vertical measurement standard deviation σ e v . Moreover, one can get a better tracking accuracy by lowering the insertion angle β or decreasing the scanning step size s.

Image processing
In the previous theoretical section, we did a complete analysis of the theoretical performance of the BC-mode OCT image. Briefly, compared with conventional B-mode OCT images, BC-mode OCT images have better SNR and better tracking accuracy. Although the image resolution deteriorates, the influence is little for samples with a flat structure. Moreover, the image resolution is on the same order of magnitude as the tracking accuracy, so it is still good enough for the surgical tool tracking task. Furthermore, the imaging speed for BC-mode OCT image visualization is fast enough for video rate surgical tool guidance.
To further improve the BC-mode OCT image quality, we will describe the image processing techniques that we use and the rationale in the following paragraphs. From the previous study, we know that for the whole signal area, the BC-mode OCT image has better SNR than the B-mode OCT image. However, signals from the surgical tool or the tissue layer, which are strong in B-mode OCT image due to the high reflectivity, become weak in BC-mode OCT image due to the averaging effect. To solve this problem, an inter frame variance processing technique is proposed to extract and enhance signals for both the surgical tool and the tissue layers.
Before proceeding, we first eliminate the saturation for all the images in the set of OCT volume data. Otherwise, the undesirable saturation line will also be extracted and enhanced. To suppress the saturation, we first decreased the optical coupling efficiency in the reference arm by adjusting the optical alignment shown in Fig. 1(a). Then we apply a Hanning window [34] to the spectrum before the Fourier transform. The Hanning window is a standard procedure for side mode suppression [35], and the process is described by Eq. (50), where I is the original spectrum,Ĩ is the new spectrum after applying the Hanning window, k refers to the k-th index in the spectrum and K is the total length of the spectrum. Next, we do the Fourier transform onĨ to get the time domain signal. The result after decreasing the optical coupling in the reference arm and applying the Hanning window to the spectrum data is shown in Fig. 9(a). We can see that the saturation is greatly reduced as well as the saturated intensity. However, there are still some saturation lines indicated by the yellow triangle in Fig. 9(a) which need to be eliminated. To detect the saturation lines, we first calculated the average intensity of each A-scan image. Then, we calculated the mean value µ and the standard deviation σ of all the average intensities. The saturation lines are identified by finding the lines whose average intensities are greater than a threshold and the threshold is calculated by T = µ + β · σ, where T is the threshold value, µ is the mean value, σ is the standard deviation, and β is a constant value which is adjustable based on how much the saturation needs to be removed. Without loss of generality, we set β = 2 in our system. After finding the saturation lines, we replace them with the lines at the same locations from the previous adjacent image slice, which is guaranteed to be free of saturation since it is processed by the same procedure before. The result after the saturation elimination is shown in Fig. 9(b). We can see that there are no obvious saturation lines.
After the saturation elimination, the main image processing scheme is shown in Fig. 10. The image set contains multiple saturation-free B-mode OCT images. The BC-mode OCT image I BC is calculated by averaging the image set, which is defined in Eq. (1). The variance image is calculated in Eq. (51), where V is the variance image, S k is the k-th image in the image set, I BC is the BC-mode OCT image, N is the total number of images in one image set and (i, j) indicates the pixel index in the image. For convenience, we rescale the variance image to [0, 255]. As we mentioned before, although the BC-mode OCT image I BC contains all the information about both the surgical tool and the tissue, the signals for both the surgical tool and the tissue layers are weakened due to the averaging effect. While on the other hand, the variance image V highlights the areas where features change fast along the Z axis, such as the surgical tool and the tissue layers. Therefore, they are good complementary to each other. To enhance the BC-mode OCT image, the idea is to add both images together. But in order to fully utilize the useful information and reduce the useless information for both images, we need to do some image processing, such as filtering, auto-threshold, and auto-rescale.Ī BC andV are the images after processing. Among them,Ī BC contains the coarse information about the surgical tool and the tissue with good contrast, whileV contains the clean detailed information such as the surgical tool and the tissue layers. The final enhanced BC-mode OCT image can be obtained by adding them together as shown in Eq. (52), whereĪ BC is the processed BC-mode OCT image,V is the processed variance image,Ĩ BC is the final enhanced BC-mode OCT image, and (i, j) indicates the pixel index in the image. The main image processing methods and results are shown in Fig. 11. For BC-mode OCT image, in order to further reduce the variance of the speckle noise, we first apply a 3 by 3 average filtering window for filtering. This step is necessary since any speckle noise will be magnified in the following contrast enhancement process, and we want the speckle noise as little as possible. Moreover, the 3 by 3 average filtering window does not affect the image resolution too much compared with the BC-mode averaging scheme. The black solid line in Fig. 11(a) show the histogram of the BC-mode OCT image after the filtering. From the histogram, we can see that the contrast of the image is not good enough since the intensities are mostly distributed from 0 to 50. This is also shown in the original BC-mode OCT image in Fig. 11(b). To further reduce the noise and enhance the image contrast, two critical intensities are found through the histogram. The first critical intensity is found at the peak of the histogram. We regard this intensity as the boundary between the noise and the signal. The second critical intensity is found through the auto triangle threshold [36]. The triangle threshold intensity is found at the point on the histogram with the largest distance to the line connecting the peak and the end of the histogram, which is indicated by the blue dotted lines in Fig. 11(a). We regard this intensity as the boundary between the coarse structure information and the detail structure information. Two critical intensities are shown as two vertical red dashed lines in Fig. 11(a), which divide the histogram into three parts. From left to right, these three parts are regarded as the noise, the coarse structure information, and the detail structure information respectively. Since the coarse structure information is the most important part in the BC-mode OCT image, we do a histogram mapping as indicated by the red dashed line shown in Fig. 11(a). The processed BC-mode OCT image resultĪ BC is shown in Fig. 11(c). Compared with the original BC-mode OCT image in Fig. 11(b), we can see that it shows the coarse structure of the cornea well with the low noise level and high contrast. However, as we discussed before, the signals of the needle and the cornea layers are weakened due to the averaging effects. But it can be solved by using the variance image, which is discussed in detail in the next paragraph.
While the BC-mode OCT image is used to obtain the coarse structure information, the variance image is used for obtaining the clean detailed structure information, such as the needle and the cornea layers. We do not use the BC-mode OCT image for that purpose since the detail information in the BC-mode OCT image is not clean enough. To reduce the noise, we first apply a 3 by 3 average filtering window to filter the variance image. The histogram of the variance image after filtering is shown as the black solid line in Fig. 11(d). To further obtain the detail information, two critical intensities are also found through the histogram. The first critical intensity is found through the auto triangle threshold [36] as indicated by the blue dotted lines in Fig. 11(d), which is explained before. We regard this intensity as the boundary between the coarse structure information and the detail structure information. The second critical intensity is found at a certain distance away from the right side of the first critical. For convenience, we set this distance the same as that between the peak intensity and the first critical intensity. In Fig. 11(d), two critical intensities are shown as two vertical red dashed lines and the histogram mapping is done as indicated by the red dashed line. The processed variance image resultV is shown in Fig. 11(e). We can see that the processed variance image clearly shows the needle and the cornea layers.
The final BC-mode OCT imageĨ BC is obtained by adding up the processed BC-mode OCT imageĪ BC and the processed variance imageV, which is shown in Fig. 11(f). The intensities over 255 are set to 255 to prevent overflow. Comparing Fig. 11(c) with Fig. 11(f), we can see that the details of the needle and the cornea layers are greatly enhanced in Fig. 11(f), while the coarse structures remain the same. All the algorithms were implemented in GPU for real-time data processing. The image and video results for the surgical tool guidance by using BC-mode OCT image visualization during microsurgery are shown in the following section.

Results
We tested the feasibility of the proposed method by guiding the insertion of a 30-gauge needle into the cornea of an ex-vivo human eye. The results are shown in Fig. 12. The three on the left are conventional B-mode OCT images. The three on the right are BC-mode OCT images. The left and the right show the same three stages during the needle insertion. The top two images show the needle just before it was inserted into the cornea surface. Middle two images show the needle when it just entered the cornea. Bottom two images show the needle when it was deep inside the cornea. We can see that in conventional B-mode OCT images, parts of the needle are missing and there are shadows underneath the needle. While in BC-mode OCT images, the whole part of the needle and the whole structure of the cornea are shown clearly. Moreover, the details of the needle, the cornea surface and the cornea bottom are well enhanced. Furthermore, it has a lower noise level with no saturation corruption. The video results of guiding a 30-gauge needle into the ex-vivo human cornea as well as pulling the needle out from the cornea using the freehand procedure can be viewed in Visualization 1. The fading of the surgical tool in Visualization 1 is mainly due to the high-frequency tremor during the needle insertion by free hand. If the surgical tool is moving at a relatively stable motion, the image of the surgical tool does not fade since the scanning will always capture the surgical tool. From both the image results and the video results, we can see that the BC-mode OCT image is beneficial for surgeons when guiding the surgical tool. For surgical procedures, surgeons can put the surgical tool roughly parallel to the densely scanning axis X at the very beginning of an operation to make it easier for the surgical tool to stay in the scanning range during the insertion. Moreover, surgeons can know the location and orientation of the surgical tool through the changes in the BC-mode OCT image when they adjust the location or orientation of the surgical tool. For example, if they try to translate the surgical tool along the sparsely scanning axis Z to see when the signal disappears, they can roughly tell if the surgical tool is in the middle or edge of the scanning area. The orientation of the surgical tool can be roughly approximated by looking at the continuity of the surgical tool signal. If the surgical tool has a large angle with the densely scanning axis X, each B-mode OCT image will capture a part of the surgical tool and this will cause a certain discontinuity of the final surgical tool signal.

Conclusions
In this paper, we propose a BC-mode OCT image visualization method for the intraoperative OCT image-guided microsurgery. Compared with B-mode OCT image, this visualization method solves both problems of shadow effects in the tissue and the difficulties of keeping the surgical tool inside the imaging field. Compared with the C-mode OCT image, this visualization method provides a faster imaging speed and larger field of view. Moreover, it provides significantly more information about the underside of the tissue surface because it shows a clearer cross-section view.
The complete theoretical analysis is developed for the performance of the BC-mode OCT image visualization, including the resolution, SNR, imaging speed and needle tracking accuracy. It is quantitatively proved that BC-mode OCT images have higher SNR and better tracking accuracy than conventional B-mode OCT images. Moreover, it has reasonable image resolution and fast imaging speed that can be used for surgical tool guidance. However, signals from the surgical tool and the tissue layers deteriorate due to the averaging effect in BC-mode OCT images. An inter frame variance processing technique is proposed to extract and enhance both signals. The enhanced BC-mode OCT image shows better signals of the surgical tool and tissue layers.
The experimental evaluation is conducted when guiding a 30-gauge needle into the cornea of the ex-vivo human eye by using BC-mode OCT image visualization. The results demonstrate that the BC-mode OCT image shows both the surgical tool and the tissue structure clearly without shadow effects. Moreover, it is easier to keep the surgical tool inside the imaging field when operating freehand by using BC-mode OCT image visualization.
Furthermore, by changing the scanning configuration, BC-mode OCT image visualization method can be applied to many other intraoperative OCT-guided microsurgery procedures. The scanning step s can be adjusted according to the tool thickness t 0 . The averaging frame number N can be set based on the necessary scanning range t and the scanning step s. However, there is a limitation to the BC-mode OCT image visualization method. When the tissue has a large curvature and the tool has a large thickness, the image resolution will deteriorate, and thus BC-mode OCT image visualization does not work well. However, for most microsurgery procedures where the tissue has a flat structure and the tool thickness is small, BC-mode OCT image visualization will be beneficial for the microsurgery guidance.

Funding
Wallace H. Coulter Foundation; Johns Hopkins University (Discovery Grant).