PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Opt Lett. Author manuscript; available in PMC 2010 June 29.
Published in final edited form as:
PMCID: PMC2893546
NIHMSID: NIHMS199024

Multiple-cardiac-cycle noise reduction in dynamic optical coherence tomography of the embryonic heart and vasculature

Abstract

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 [5].

Several techniques for reduction of coherent and incoherent noise in OCT images have been used [6], including cross-correlation image alignment followed by weighted averaging [7], compounding several B-scans acquired at neighboring positions in the sample [8], and elastic spatial registration of several OCT B-scans followed by averaging [9]. Transform-domain methods, which are based on the assumption that the signal (but not the noise) can be represented sparsely in wavelet [10] or curvelet [11] function families, have been applied successfully to piecewise smooth retinal OCT images. Hardware-based methods have also been reported [12]. 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 [2]. 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 [11]). Accordingly, we model the measured images as

Im(x,t)=Ip(x,t)+It(x,t),
(1)

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

Ip(x,t)Ip(x,t+T),
(2)

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<tT, 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

Ic(x,r,t)=Im(x,T+(r1)σT+t),0t<L,
(3)

with σT>LTmin=(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).

Fig. 1
Multicycle noise reduction: (a) raw sequence, (b) cutting, (c) temporal registration, (d) noise reduction via sample median along the r dimension.

Next, the cut sequences Ic are warped so as to match the template sequence I0. We define the cost functional Q:MR+ as

Q{wr}=(1λ)0TI0(x,t)Ic(x,r,wr(t))dtdx+λ0Tddtwr(t)1dt.
(4)

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 M={wC1([0,T))0w(t)<Landw(t1)<w(t2),t1<t2} 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 [13]. The optimal warping function produced by this algorithm is wr{wrMQ{wr}=minwMQ{w}}.

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,wr(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,

I~p(x,t)=Median({I0(x,t),Ic(x,r,wr(t))1rR}),0t<T.
(5)

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

Fig. 2
(Color online) Frames from SS-OCT of (a) raw, (b) 3×3 spatial moving average, (c) four-frame temporally binned, and (d) noise-reduced (via proposed multicycle sample median) sequences of the beating heart in an E9.5 mouse embryo. Arrow indicates ...

We quantified the results presented in Fig. 2 using several image-quality metrics initially proposed to evaluate OCT noise-reduction algorithms on 2D images [10] that we adapted to our 2D+time data: (i) the global signal-to-noise ratio, SNR=10log[max(Xlin2)σ12], (ii) the contrast-to-noise ratio, CNR=1(M1)m=2M10log[(μmμ1)(σm2+σ12)12], and (iii) the equivalent number of looks, ENL=1(K1)m=2K(μm2)σm2, 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 [10], β, 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 β.

Table 1
SNR, CNR, and ENL of Data in Fig. 2

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 [14]. 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.

Acknowledgments

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

Footnotes

OCIS codes: 110.4500, 180.1655, 030.4280, 030.6140.

References

1. Fleming C, Ripplinger C, Webb B, Efimov I, Rollins A. J. Biomed. Opt. 2008;13:030505. [PMC free article] [PubMed]
2. Larina I, Sudheendran N, Ghosn M, Jiang J, Cable A, Larin K, Dickinson M. J. Biomed. Opt. 2008;13:060506. [PubMed]
3. Larina I, Ivers S, Syed S, Dickinson M, Larin K. Opt. Lett. 2009;34:986. [PMC free article] [PubMed]
4. Liu A, Wang R, Thornburg K, Rugonyi S. J. Biomed. Opt. 2009;14:044020-1–11. [PubMed]
5. Gargesha M, Jenkins MW, Rollins AM, Wilson DL. Opt. Express. 2008;16:12313. [PMC free article] [PubMed]
6. Schmitt J, Xiang S, Yung K. J. Biomed. Opt. 1999;4:95. [PubMed]
7. Sander B, Larsen M, Thrane L, Hougaard J, Jørgensen TM. Br. J. Ophthamol. 2005;89:207. [PMC free article] [PubMed]
8. Popescu DP, Hewko MD, Sowa MG. Opt. Commun. 2007;269:247.
9. Jørgensen TM, Thomadsen J, Christensen U, Soliman W, Sander B. J. Biomed. Opt. 2007;12:041208. [PubMed]
10. Adler DC, Ko TH, Fujimoto JG. Opt. Lett. 2004;29:2878. [PubMed]
11. Jian Z, Yu Z, Yu L, Rao B, Chen Z, Tromberg B. Opt. Lett. 2009;34:1516. [PMC free article] [PubMed]
12. Wang H, Rollins A. J. Biomed. Opt. 2009;14:030512. [PMC free article] [PubMed]
13. Liebling M, Vermot J, Forouhar A, Gharib M, Dickinson M, Fraser S. IEEE International Symposium on Biomedical Imaging; 2006.p. 1156.
14. Braun S. Acustica. 1975;32:69.