|Home | About | Journals | Submit | Contact Us | Français|
Recent progress in optical coherence tomography (OCT) allows imaging dynamic structures and fluid flow within scattering tissue, such as the beating heart and blood flow in mouse embryos. Accurate representation and analysis of these dynamic behaviors require reducing the noise of the acquired data. Although noise can be reduced by averaging multiple neighboring pixels in space or time, such operations reduce the effective spatial or temporal resolution that can be achieved. We have developed a computational postprocessing technique to restore image sequences of cyclically moving structures that preserves frame rate and spatial resolution. The signal-to-noise ratio (SNR) is improved by combining images from multiple cycles that have been synchronized with a temporally elastic registration procedure. Here we show how this technique can be applied to OCT images of the circulatory system in cultured mouse embryos. Our technique significantly improves the SNR while preserving temporal and spatial resolution.
Fast optical coherence tomography (OCT) is being used increasingly for imaging of embryonic cardiovascular dynamics in animal models, as it is noninvasive and allows imaging through several millimeters of biological tissue with single-cell resolution [1-4]. To characterize the highly dynamic cardiovascular system, it is critically important to develop effective noise-reduction methods that preserve spatial and temporal resolution without introducing blurring .
Several techniques for reduction of coherent and incoherent noise in OCT images have been used , including cross-correlation image alignment followed by weighted averaging , compounding several B-scans acquired at neighboring positions in the sample , and elastic spatial registration of several OCT B-scans followed by averaging . Transform-domain methods, which are based on the assumption that the signal (but not the noise) can be represented sparsely in wavelet  or curvelet  function families, have been applied successfully to piecewise smooth retinal OCT images. Hardware-based methods have also been reported . Drawbacks of the above techniques include loss of spatial or temporal resolution, reliability on object sparsity and geometry assumptions that may not apply to all image types, or modification of the instrumentation itself.
Here, we propose a technique that takes advantage of the periodic (or nearly periodic) organization of cardiovascular image sequences along the temporal axis, which provides a natural framework for noise reduction in which neither spatial smoothness, object geometry assumptions, nor optical setup modifications are required.
We consider an input image sequence of a cardiovascular structure (e.g., the heart or the aorta) whose grayscale values are proportional to the amount of scattering in the tissue or blood as measured by a swept-source OCT system . We assume that the imaged structures either are static or undergo a periodic motion that follows the heartbeat. We also assume that additive noise contributing to the image at each space–time location follows a zero-mean random distribution and is uncorrelated with the signal (to fit this assumption in the case of multiplicative noise, the raw data can be processed with a logarithmic transform prior to the proposed noise-removal method and postprocessed with an exponential transform ). Accordingly, we model the measured images as
where x=(x,y), Ip(x,t) is the contribution from the permanent structures (such as the heart wall), and It(x,t) is the intensity contribution from the transient structures (such as additive noise or spurious scatterers that are neither static nor moving cyclically). Our assumption of permanent structures contributing to the image contrast in a periodic fashion translates to
where T is the period of the cardiac cycle.
Our algorithm aims at estimating Ĩp(x,t), the image contribution from the periodic features, given the measured intensity Im(x,t). We require Im to be measured over multiple periods, that is, 0≤t<NT, with N the total number of periods. We consider the first heartbeat period, I0(x,t)=Im(x,t) with 0<t≤T, as a template to be matched to subsequent periods. We assume that the length of the latter periods may vary by an estimated percentage α such that the minimal and maximal period lengths are Tmin=(1–α)T and Tmax=(1+α)T, respectively. We next extract R shorter sequences from the original measured sequence. To make sure that it contains at least one complete cycle of the heartbeat that is in phase with the template, each extracted sequence must be of length L≥2Tmax. However, since the template is matched over only a portion of maximal length Tmax, the rest of the data is wasted if the extracted sequences are disjoint. We therefore allow for some overlap when extracting sequences, yet we must ensure that the template is not matched to the same region in two different extracted sequences. This is guaranteed as long as the overlap does not exceed Tmin. Given the above conditions, we define the extracted sequences as
with σT>L−Tmin=(1+3α)T as the step size controlling the overlap of the cut sequences, and r=1, … ,R the cut index, where R is such that T+(R−1)σT+L<NT (see Fig. 1).
Next, the cut sequences Ic are warped so as to match the template sequence I0. We define the cost functional as
The first integral compares a temporally warped sequence (a sequence whose temporal axis is deformed) with the template, whereas the second integral keeps the extent of this deformation in check by penalizing warping functions that stretch or compress the time axis. These two contributions to the cost function are balanced by the parameter 0≤λ<1 to favor either good matching of the warped and reference sequences (λ=0) or the temporal integrity of the warped sequence (λ→1). We note that the warping functions wr are in the set of continuous, nonnegative, and monotonically increasing functions bounded by L and defined over the interval [0, T). We then determine a warping function that minimizes the cost function. This minimization is carried out using a previously described dynamic programming algorithm . The optimal warping function produced by this algorithm is .
In the registered sequences, signal contribution from permanent (static or cyclically moving) structures occurs at the same space–time position in every sequence. Corresponding pixels in each of the R sequences, Ic(x,r,r(t)), therefore represent as many measurements of the same spatial configuration, with an additive random noise contribution. We can therefore estimate the contribution from permanent structures Ip to the measured signal Im by finding the sample median along r,
Importantly, this operation reduces the noise floor as it virtually increases the integration time while preserving structural information. Choosing the sample median rather than an average ensures robustness to outliers, that is, sparse contributions from any nonperiodic signal.
To verify our method in practice, we acquired image sequences of a live 9.5 day mouse embryo using a swept-source OCT (SS-OCT) system as described previously in . Figure 2(a) shows a representative frame of an SS-OCT image series of the embryonic heart and yolk sac acquired at 50 frames per second (fps), each of size 256×512 pixels (512 pixels per A-scan). A 3×3 moving average filter applied spatially reduces the graininess of the images but at the expense of spatial resolution [Fig. 2(b)]. Averaging four successive A-scans reduces temporal resolution [Fig. 2(c)]. Figure 2(d) shows the result of our proposed multicycle method, based on the 1000 frameslong sequence of Fig. 2(a) containing approximately 40 heartbeats and resulting in R=26 shorter sequences over which the sample median was determined. We assumed T=0.5 s and α =15% and set σ=1.5, L=2.5Tmax, and λ=0.5. Our technique reduces the noise floor without perceptible loss in spatial or temporal resolution. Interestingly, speckle is preserved in some locations in Fig. 2(e). This can be explained by the absence of nonperiodically moving structures in the optical path leading to these locations that would make the speckle pattern temporally aperiodic.
We quantified the results presented in Fig. 2 using several image-quality metrics initially proposed to evaluate OCT noise-reduction algorithms on 2D images  that we adapted to our 2D+time data: (i) the global signal-to-noise ratio, , (ii) the contrast-to-noise ratio, , and (iii) the equivalent number of looks, , where Xlin is the 3D matrix of pixel values in the OCT image sequence on a linear intensity scale, μ1 and σ1 are pixel mean and standard deviation (SD) of a background region in the image sequence, μm and σm are the pixel mean and SD of the mth 3D ROI in the image sequence. CNR and ENL values are computed after logarithmic transformation of the image sequence. Three homogeneous ROIs, m=2, … ,K, with K=4 are defined in regions without edges for computing the ENL. The temporal extents of these ROIs are defined so that the average pixel intensity variation in these regions remains constant. Seven ROIs (homogeneous and containing edges; m=2, … ,M, M=8) are used for the CNR. The first ROI, m=1, is used to compute the background-noise level. The 3D ROIs are shown in Fig. 2(e) and 2(f). In addition, a global image sharpness parameter , β, is used to determine how much the noise-reduction process has degraded the edge sharpness. A β value close to unity means edges are preserved. The metrics computed for the sequences in Figs.2(a)–2(d) are shown in Table 1. Our proposed technique produces significant improvement of SNR, ENL, and CNR, confirming the visual observation that our technique improves the SNR while preserving edges, as indicated by β.
In summary, we have validated a noise-reduction algorithm for OCT images of cardiovascular structures whose motion is cyclical, as dictated by the beat of the heart. The essence of this technique, combining images from multiple cycles that have been phase locked to obtain an estimate of the noise-free signal, is reminiscent of methods employed for extracting strictly periodic monodimensional waveforms . Similarly, the SNR can be increased as the number N of imaged heartbeat cycles is increased, as long as the motion of the structure is stereotypical. Here, in order to take into account slight variations in the heart rhythm, registration is nonrigid along the temporal dimension (yet no spatial registration is carried out to preserve the original structure shape). Our technique does not require assumptions on the spatial smoothness, organization, or shape of the structures being imaged and therefore allows preserving both temporal and spatial resolution.
This research was supported by Department of Defense Office of Naval Research Young Investigator Program (K. V. L.), National Institutes of Health (NIH) HL077187, HL095586, and T32HL007676, and a UCSB Regents Junior Faculty Fellowship (M. L.).
OCIS codes: 110.4500, 180.1655, 030.4280, 030.6140.