Lateral resolution improvement of oversampled OCT images using Capon estimation of weighted subvolume contribution

: A novel technique for lateral resolution improvement in optical coherence tomography (OCT) is presented. The proposed method is based on lateral oversampling of the image. The locations and weights of multiple high spatial resolution sub-volumes are calculated using a Capon estimator assuming each contributes a weighted portion to the detected signal. This technique is independent of the delivery optics and the depth of field. Experimental results demonstrate that it is possible to achieve ~4x lateral resolution improvement which can be diagnostically valuable, especially in cases where the delivery optics are constrained to low numerical aperture (NA).


Introduction
The lateral resolution of Optical Coherence Tomography (OCT) images is fundamentally limited by the delivery optics. Although a higher effective numerical aperture (NA) can provide better lateral resolution, it limits the depth of field (DOF) of the imaging. The mutually exclusive relationship between lateral resolution and DOF imposes significant limitations and, thus, design trade-offs become necessary. These limitations become even more pronounced when the sample beam focusing is physically constrained, e.g. in the case of ophthalmic OCT imaging where the optics of the eye itself define the resolution. Furthermore, the lateral resolution degrades away from the focal plane blurring the images when moving out-of-focus.
A number of different techniques have been developed in order to improve the lateral resolution or extend the depth of imaging without resolution degradation. Adaptive optics [1] or axicon lenses [2] have been used to improve or maintain the resolution over larger scanning depths. Dynamic focusing [3] and focus tracking [4] in the sample arm have also been employed. However, in all of these cases special hardware is required which limits the scanning speed and the applicability to real-time or in vivo imaging. The lateral resolution away from the focus can also be numerically corrected using inverse scattering algorithms [5,6]. Interferometric synthetic aperture microscopy (ISAM) can achieve depth-independent resolution throughout a volume even when the focus is fixed at one depth [7]. The disadvantage of this method is a decrease in signal-to-noise ratio (SNR) away from the focus.
A novel two-dimensional numerical method to achieve improved lateral resolution and a long measurement range was also developed [8]. However, this method is extremely sensitive to the phase stability of the measurements, which may be easily maintained during in vivo imaging where sample movement is unavoidable. Deconvolution methods have also been used for lateral resolution improvement, usually limited to a factor of < 2x enhancement. A noniterative numerical method for laterally superresolving OCT images, using a depth dependent PSF [9] and Gaussian beam deconvolution were also described [10,11]. Wiener and Lucy-Richardson OCT image restoration algorithms were implemented and compared proving that the later provides better performance [11]. Deconvolution-based methods usually require that the point spread function (PSF) of the system is known a priori. More recently, an automatic PSF estimation technique to deconvolve the OCT images, based on the discontinuity of information entropy was demonstrated [12].
In this paper, we present a new method for lateral resolution improvement which is independent of the focusing conditions and the depth of field. It is based on lateral oversampling [13,14] and the estimation of the weighted contributions of multiple signals forming high spatial resolution sub-volumes, based on Capon estimation. This approach is inspired from radar range oversampling techniques [15,16]. The PSF is also required for this technique but it can be estimated from the image using blind deconvolution. Results from the implementation of this method on laterally oversampled OCT images have shown that it is possible to achieve lateral resolution improvement of about a factor of 4, greater than that achieved with any of the deconvolution approaches.

Lateral oversampling in OCT
If the lateral resolution of an OCT system is Δx then, according to the Nyquist sampling theorem, the information in the transverse direction can be completely determined if the distance between the centers of two adjacent resolution volumes is Δx/2. If the lateral direction is oversampled by a factor of L, with a sampling rate of L/Δx, adjacent resolution volumes will overlap, since they are shifted by Δx/L, and the signals from successive volumes will be correlated. This is shown schematically in Fig. 1 for L = 3 and a square PSF of width Δx. Fig. 1. A schematic diagram of L = 3 lateral oversampling. Each oversampled signal Vi consists of independent signals from L subvolumes (S i to S L + i-1 ) with a 2/3 overlap between volumes which share signals Si to S L + i-2 with the previous volume [15]. The red rectangle indicated the shared volume which can be isolated from appropriately combining the oversampled signals.
Using the theory of wave scattering from a random medium, Zhang et al. (2005) described the possibility for resolution refinement from oversampling signals [16]. The resolution improvement stems from the fact that the most significant contribution to the crosscorrelation function, of signals from two overlapping resolution volumes, is from scatterers in the shared volume. Signals from all overlapped volumes can be combined optimally to isolate the contribution of the scatterers in that narrow common region at the center of the overlap (red rectangle in Fig. 1) thus improving the resolution without actually breaking the physical resolution limit set by the NA of the optical system.
If a resolution volume of width Δx is assumed to include L high spatial resolution subvolumes, each of a size Δx/L, then the OCT signal at each lateral location, V i (t), is a weighted contribution of those L independent high spatial resolution signals (S i to S L + i-1 ). The weight of each subvolume is determined by the PSF of the system. The relationship between the high spatial resolution and the OCT signals can be represented by the following matrix notation: is a column vector of oversampled signals, V i (t), from L successively spaced resolution volumes, [ ] is a vector of independent high spatial resolution signals (which could be individual or continuous scatterers) at 2L-1 sub-volumes, within the range covered by L consecutive resolution volumes, and the superscript T is the transpose. Moreover, [ ] is a matrix of shifted PSFs and has a size of L(2L-1). To more clearly illustrate the meaning of Eq. (1), consider the simplified case of Fig. 1 (L = 3 oversampling and a square PSF). In this case, Eq. (1) takes the following form: [ ] Each oversampled signal is simply a summation of equally weighted independent signals from those sub-volumes contained in each resolution volume.
Equation (1) The above formulation depends and will yield the signal from the high-resolution subvolumes, Si, which may contain a single or multiple scatterers or even a continuous interface. The resulting Si will be the combined effect of the contents of the sub-volume. This is not unlike the signal that would be obtained using a narrower beam to image that same subvolume independently. However, since, in the proposed approach, the adjacent sub-volumes are simultaneously illuminated there is a chance for interferometric interference which would not be present in the case of imaging with a narrower beam. Fortunately, since OCT mainly detects backscattered signals, the chance of interference between the reflection of laterally displaced scatterers is minor, as also shown by similar models in the literature [5,6,15,16]

Resolution enhancement using an inverse matrix solution
In order to reconstruct an OCT image with improved lateral resolution, the high spatial resolution sub-volumes Si must be determined. One approach is to extract those sub-volumes by directly solving Eq. (1), where A −1 is the inverse of matrix A. However, the divisions involved in the formation of the inverse matrix, make this approach extremely sensitive to noise (noise amplification) and limit the practical applicability of the technique. This will become immediately apparent when the experimental results are described in the next section.

Resolution enhancement using Capon's method
A more robust approach to solve the inverse problem and recover the high spatial resolution signals, S i , is to find the optimum weights which relate the V i 's and the S i 's. Capon's method is one approach which can be used to find the appropriate weight function. This technique was extensively studied in the literature especially for high-resolution spectral estimation, including Doppler ultrasound [17][18][19]. Capon's approach to spectral estimation is described first, followed by its adaptation for OCT lateral resolution enhancement.

Capon spectral estimation
The Capon spectral estimation method is a filter-bank approach where the problem of filter design is solved with some specific constraints. These constraints are that the signal at a specific frequency, ω, is passed through unchanged (unit gain) and that the output power of the overall frequency domain is minimized. Appropriate choices for these constraints ensure that the Capon filter is actually a matched filter [20,21]. Let be a finite impulse response (FIR) filter of M-th order, where H denotes conjugate transpose. If x(n) is the input data sequence, then the output of the filter at time n is given by: If R is the covariance matrix of the data vector ( ) x n such that where E is the expectation operator, then the power of the filter output is The filter frequency response is given by where a(ω) is the steering vector defined as: and T denotes the transpose. Capon's method uses a bandpass filter satisfying: subject to Mathematically this is equivalent to minimizing the following cost function: where μ is a Lagrange multiplier. The minimization of this equation leads to the solution:

Adaptation of Capon's method to laterally oversampled OCT signals
The oversampled OCT signal, including noise, can be expressed as: is a matrix of weighted filters with a size of L(2L-1), and each of its columns is the corresponding weighting vector, and is the output of the weighted filter, corresponding to an estimate of the high spatial resolution sub-volume signal, the following expression is obtained: Obviously, if the initial high resolution signal Si could be perfectly reconstructed, w i a i should be equal to unity while the other terms should be minimum. Using Capon's approach, the weights can be obtained by solving a constrained optimization: The power to be estimated at sub-volume i, is a weighted sum of the power from all 2L-1 subvolumes. The optimization ensures that the contamination, from sub-volumes other than the i th sub-volume is minimized while the response of the weighting vector at the sub-volume in which the estimation is performed is unity. Using a Lagrange cost function, the weighting vectors can be obtained as follows: is the autocorrelation matrix of the oversampled signals and each element is: which is the cross correlation function of the i th and j th oversampled signal and a i is the PSF for i. The process of reconstruction is shown graphically in Fig. 2 where five high resolution volumes (S 1 -S 5 ) are reconstructed from three oversampled signals (V 1 -V 3 ) using the weights (w 11 -w 35 ) estimated using Capon's approach. Fig. 2. Schematic illustrating graphically the process of lateral resolution enhancement using oversampling and Capon's method to reconstruct the high spatial resolution signal. In this example, five high resolution volumes (S 1 -S 5 ) are reconstructed from three oversampled signals (V 1 -V 3 ) using the weights (w 11 -w 35 ) estimated using Capon's approach.

PSF estimation
The PSF is necessary to extract the Capon weights. One approach could be to directly measure the PSF from an infinitesimally small object. However, this is experimentally challenging. Alternatively, the PSF can be extracted from the image using optimization techniques. In this work, the PSF of the system was estimated using blind deconvolution based on the iterative Richardson-Lucy algorithm [22]. Blind deconvolution can be used effectively with no a priori information about the distortion (blurring and noise) of the image.
where O n + 1 and O n are the current and previous estimate of the desired image, P is the PSF, P n + 1 and P n are the current and previous estimate of the PSF and P and n O are the complex conjugate of P and O n . I is the original blurred image. The algorithm begins with an estimate for both P and O, and, by constraining those, tries to obtain the original image. The constrains of the algorithm are positivity and stability for both the image and the PSF.
An important characteristic of this method is that when the SNR of the original image is low, the PSF estimate is not as accurate. The implication, for the purposes of resolution improvement with either deconvolution or the proposed Capon-based method, is that the resolution improvement is proportional to the SNR of the original image.

Results and discussion
The experimental evaluation of the methods described above was performed on OCT images from a phantom of microspheres with an oversampling factor of 28. A swept-source system operating at 1.3 μm center wavelength with a lateral resolution of 14μm was used to acquire 2000 A-scans per image at a step of 0.5μm. Furthermore, images from skin, with the same oversampling factor, were also acquired and processed to verify the applicability of the proposed method on biological samples with noise and speckle. In addition to the proposed Capon-based method, the inverse matrix solution and deconvolution were also performed for comparison purposes.
The PSF of the system was estimated using blind deconvolution at one depth of the image and the same PSF was then applied to the whole image. This approximation was considered adequate for the case of the low NA imaging system that was used in this study. Figure 3 illustrates the validity of this assumption since PSFs estimated at different depths did not differ by more than 5 μm. However, for high NA imaging, the PSF should be estimated separately for each depth in the sample. In order to estimate the lateral resolution of the images, the full-width-at-half-maximum (FWHM) of one peak, which corresponded to a single microsphere, was measured. However, since the intersection of the scattering particle with the imaging volume would not, necessarily, be through the center of the sphere, the diameter of the scatterer could not be expected to have the original sphere diameter of 6 μm. Instead, the system resolution was measured, using the knife edge method, and that width was subtracted from the FWHM of the peak to estimate its effective diameter. This approach would also be applicable to spheres larger than the PSF. The result is, of course, an approximation since the cross-section of a sphere partially covered by the imaging beam is no longer a sphere and the PSF after convolution with the scatterer may not result in a peak width of exactly the sum of the two (depending on the shape). However, this is an adequate for first order comparisons of the resolution improvement from the various approaches.

Inverse matrix solution
Images from a phantom of 6 μm microspheres embedded in acrylamide gel, were acquired with an oversampling factor of 28. Subsequently they were processed as described above. Matrix A was formed using the estimated PSF at a depth ~300μm below the surface of the sample, and the high resolution image was obtained by direct solution of Eq. (5). The results of the inverse matrix solution are shown in Table 1 and Fig. 4. The comparison of the images confirms that the resolution after processing improved. However, there was significant noise amplification. In order to reduce the noise, filtering was applied to the images, which in turn degraded the resolution. For standard OCT, the FWHM of the single sphere was 17.85 μm, for the OCT image after deconvolution, 12.85 μm, and for the inverse matrix solution with filtering, the FWHM of the same peak became 9.38 μm. The resolution could have been finer if less filtering was performed on the final image but, in that case, there would have been significantly more noise in the image. The FWHM of the signal from a single peak is determined by the lateral PSF of the system plus the effective width of the sphere. For quantification purposes, considering that the lateral resolution of the system is 14 μm then the width of the sphere, at the location of measurement of the FWHM, is 17.85-14 = 3.85μm. By subtracting from each measurment the effective width of the sphere, the PSF of the system was estimated and the resolution improvement quantified. where I max is the maximum value of the image, and σ b is the standard deviation of the noise in the from a portion of the image where there is only background noise. The resolution improved by a factor of ~2.5 with the inverse solution method. This is a slightly better result than deconvolution which improved the resolution by a factor of 1.55. The main limitations of the inverse solution method are that it is very sensitive to noise and that the performance of this approach degrades when the original image has a low SNR. The computational time for each image was ~15.6 seconds, on a 2.5GHz PC using Matlab, for the inverse matrix solution and ~6.2 seconds for the deconvolution of the OCT image. For the inverse solution, the most computationally intensive step was the calculation of the inverse matrix.

Capon weight method
The Capon weight method was applied to the same oversampled images as before. This method was implemented using three different approaches each time employing a cruder approximation of the weights in order to reduce the computational time.

Capon weight method with localized weight calculation
In the first approach, for each line at a time and repeating for all lines, the value of each subresolution volume (pixel) was estimated based on a weight vector calculated from the values of a 56 pixel (28 μm) neighborhood around the center pixel. An example is shown in Fig. 5. Table 2 and Fig. 6 provide illustrate, visually and quantitatively, that the Capon weight method provides excellent improvement with no noise amplification after processing. Using this approach, the resolution improved by a factor of 3.73. This is a significant improvement, compared to deconvolution where the resolution improved by a factor of 1.55. The SNR, after applying the weighting method, did not decrease since the noise was not amplified. This method is computationally complex. It took ~1363 seconds to process an image on a 2.5GHz PC. This was approximately 100 times longer than the inverse matrix solution. The most computationally intensive step here was the calculation of several inverse matrices for each line of the image.

Capon weight method with per line weight calculation
In the second approach, for each line at a time and for all lines, the weight vectors for all subvolumes were calculated from the values of all pixels of the line, and used to estimate the values of each sub-resolution volume (pixel) individually. The results are shown in the Fig. 7 below. Table 3 summarizes the results. In this case, the resolution improved by a factor of 3.37. This was, again, superior compared to deconvolution. However, the SNR after the Capon weighting decreased in this case. In addition, this approach was the most computationally complex. It took ~8023 seconds to process an image on a 2.5GHz PC. This time was approximately 500 times longer than that of the inverse matrix solution. The most computationally intensive step here was the calculation of an inverse matrix of 2000x2000 data points for each line of the image.

Capon weight method with single line weight estimate
In the third Capon-based approach, the weights of all pixels for one line were computed as in the previous section but the same weights are used for all lines. The results are shown in Table 4 and Fig. 8 below. It is clear that there is significant resolution improvement (2.82x) compared to deconvolution (1.55x). In addition, the SNR was not affected in this case. This approach was less computationally intensive than the previous two. It required ~30.4 seconds to process an image on a 2.5GHz PC, only approximately 2 times longer than that for the inverse solution. The most computationally intensive step here was also the calculation of an inverse matrix but it was calculated only one time per image. The algorithm was also tested on human tissue images to evaluate the performance of the method on data with speckle and noise. The results of imaging normal and invasive adenocarcinoma of the colon during surgical excision are shown in Fig. 9 below. The resolution improvement is evident when examining the gland and crypts structures in the zoomed portions of the images. It also obvious that speckle does not affect the performance of the proposed technique.

Conclusions
A novel technique for improving the lateral resolution of oversampled OCT images was developed and experimentally verified. The improvement is independent of the focusing optics and the depth of field. The new technique is based on calculating high spatial resolution sub-volumes from overlapping OCT A-Scans. Unlike the inverse matrix solution of the problem, using the proposed Capon weight method, the resolution can be improved without noise amplification and in the presence of speckle. Applying the weighting algorithm to overlapped and shifted by one sub-volume regions of the signal, the resolution improved by a factor of ~4, and the SNR was not affected. The computational load of the technique can be significantly reduced by computing the weights only one time and applying the same weights to all lines. However, in this case, the resolution improves only by a factor of ~3. It should be noted that deconvolution can only improve the resolution by < 2. The performance of the proposed technique could decrease for low SNR images since it depends on an accurate estimate of the PSF which in turn depends on the SNR. In addition to the increased computational complexity, this method also increases the acquisition time since it requires oversampling in the lateral direction. However, given the progress in current OCT technology, both in terms of imaging speed and GPU-based computational implementations, application of the proposed technique for real-time imaging is feasible and could provide considerable diagnostic advantages as a result of the much needed resolution improvement.

Funding
This research was partially supported by the Research Promotion Foundation of Cyprus and the KIOS Research Center for Intelligent Systems and Networks.