Speckle suppression via sparse representation for wide-field imaging through turbid media

Hwanchol Jang; Changhyeong Yoon; Euiheon Chung; Wonshik Choi; Heung-No Lee

doi:10.1364/OE.22.016619

1. Introduction

Coherent waves propagating through a turbid medium experience multiple scattering and the output waves have various propagation directions and phases. Due to these variations, the outgoing waves are spatially scrambled and become a speckle field (SF) at an observation plane; a speckle refers to an interference pattern with random spots overlaying an image. The SF is incoherent to the original wave and it was regarded that the information of the original wave is lost. However, recent studies show very interesting results that the multiple scattering in the medium can be discarded, or even be utilized for higher resolution purpose in optical imaging systems [1–8].

The wave propagation is a time reversible (TR) process and the multiple scattering can be reversed by a TR operator. The inversion is done for i) image recovery [1–7] which retrieves the original input waves from the output SF, and for ii) wavefront shaping [1,8] which manipulates the input waves so that the output waves have a desired form. These are the two main branches of the current studies relating the optical imaging through turbid media [1]. Phase conjugation (PC) is the monochromatic version of the TR operator and it is used in optical imaging through turbid media for speckle suppression [1,2]. PC compensates the phase variations due to multiple scattering in turbid media by recording the SFs and back-propagating the complex conjugates of them through the media so that the phase variations are cancelled.

This PC can be done virtually through computational estimation in digital signal processing rather than through optical compensations if the so-called transmission matrix (TM) [3] whose columns are the output SFs of the medium for a set of input basis waves is available. It has been demonstrated that the computational PC through estimation in the TM-based approach can reverse the multiple scattering from the turbid media both in the image recovery and the wavefront shaping; and the use of the PC is common for speckle suppression in the TM-based approach [1,3–8]. Here, the (digital) TM-based approach has a number of advantages over the optical approaches for it has an image data format which is reproducible [4].

From now on, we restrict our discussion to the TM-based image recovery and use computational version of PC as a point of reference. For simplicity, the PC will stand for computational PC. Our research starts with an observation that PC in the TM-based recovery is very delicate in the presence of correlations among the columns of a TM. During inversion, the SF in the observation plane which is correlated to those of more than one basis waves is removed only partially, and the recovered image by PC still suffers from the remaining speckle [6,9]. Additional speckle suppression is needed and this requires ensemble averaging over many samples.

In this paper, we aim to reduce the speckle in the TM-based recovery process for high resolution imaging through a turbid medium. We find that the main cause of the remaining speckle in the TM-based recovery is the error in the estimation. We apply for a relatively new signal reconstruction framework referred to as sparse representation (SR). The SR has been reported to provide superior estimation even in ill-conditioned systems [10], and found to be very useful in holography and microscopy applications [11]; note that a special case of SR for measurement reduction is the well-known compressed sensing. We suppress the speckle in the TM-based recovery through SR. For this, we show the validity of the SR framework in the context of the TM-based recovery. We show that the speckle in the TM-based recovery can be suppressed effectively via SR.

Recently, many imaging applications such as microscopy [11], [12], holography [13], and MRI [14] deal with underdetermined systems where the dimension of the measurement is smaller than that of the signal [11]. The rationale behind this is that the information contents of most natural objects such as living organisms are rather small to be stored in the measurements whose dimension are larger than or equal to the number of elements of the representing basis, such as the Dirac basis (pixels) or the Fourier basis (angular spectral elements), in which we want to demonstrate the objects, and can be contained in much smaller dimensional measurements [11]. We thus consider an underdetermined measurement system in this paper for the same reason. Besides, it is definite that the measurement dimension reduction is appropriate for efficient resource utilization especially in the TM-based imaging which requires large amount of data acquisitions.

This paper is organized as follows: Section 2 describes the TM-based image recovery, and Section 3 presents the PC method and discusses the problem of it in the TM-based image recovery. Section 4 presents the SR framework for speckle suppression in the TM-based image recovery. Experimental results are discussed in Section 5, and Section 6 concludes the paper.

2. System description

2.1 Experimental set-up

The baseline method considered here is the turbid lens imaging (TLI) which is a digital TM-based recovery scheme [5,9]. Consider the experimental schematics in Fig. 1 [5,9]. The turbid medium is placed between the sample plane (SP) and the objective lens. It is a 25 μm thick layer of ZnO nanoparticles deposited on a slide glass. In the experiment, a collection of plane waves ( $λ = 633 nm$ ) each with different incident angle to the turbid medium is used as a basis ([Fig. 1(a)]); the incident angle to the turbid medium is controlled by the galvanometer. The SF for each plane wave is obtained and stored as a column in TM ([Fig. 1(b)]). Thus, TM is a matrix whose columns are the SFs. For an object, we write patterns similar to USAF target on a spatial light modulator (SLM). The object SF (OSF), which is the output SF of turbid medium with the object wave, instead of the plane wave, is then obtained ([Fig. 1(c)] and [Fig. 1(d)]). Note that all the recorded images are processed as E-field images by Hilbert transformation [15]. With the TM and the OSF, the estimation in [5,9] is done by PC digitally. There are other TM-based approaches using PC worth to make note of, such as the work by Popoff et al. [3], which uses different basis waves from ours, and the work by Hillman et al. [7], which employs a different experimental configuration.

Fig. 1 Experimental schematics. (a) Recording of the TM. SP: sample plane, BS: beam splitter, SB: sample beam, RB: reference beam. (b) The recorded TM. Only the intensities of the SFs are shown here. (c) Recording of the OSF. SLM: spatial light modulator. (d) The recorded OSF. Only the intensities of the SFs are shown here. Scale bar: 10 μm.

Download Full Size | PDF

2.2 Image recovery using transmission matrices

Consider an optical wave containing an object image. It can be decomposed into a set of plane waves with various angles of propagation, formally written as

o (x, y) = \sum_{k_{x}} \sum_{k_{y}} a (k_{x}, k_{y}) p (x, y; k_{x}, k_{y})

where

o (x, y)

is the object wave,

a (k_{x}, k_{y})

and

p (x, y; k_{x}, k_{y})

are the angular spectrum and the plane wave, respectively, of the object wave with the propagation angle with k_x and k_y; k_x and k_y are the wave vector components of the wave in the x and y directions. Each plane wave with a distinct incident angle to the turbid medium has a unique output SF at the observation plane. As a result, the OSF is expressed as

y (x, y) = \sum_{k_{x}} \sum_{k_{y}} a (k_{x}, k_{y}) t (x, y; k_{x}, k_{y}) + n (x, y)

where

t (x, y; k_{x}, k_{y})

is the output SF of the medium with the input plane wave

p (x, y; k_{x}, k_{y})

and

n (x, y)

is the additive noise. With

t (x, y; k_{x}, k_{y})

for all

(k_{x}, k_{y})

considered, the estimate of the angular spectrum,

\hat{a} (k_{x}, k_{y})

, is made for a given

y (x, y)

. Using the angular spectrum estimation, the recovery of the object image is made. Severe speckle is observed in the recovered object image [6,9]. The speckle becomes the major obstacle for object image interpretation. The common method to mitigate the speckle is to generate many object samples with uncorrelated speckles and to take the average of them [16]. This method is most commonly used in the TM-based recovery [6,9]. But, this averaging lengthens the image acquisition time and thus is not desirable. In time-critical cases or in the case of imaging a moving object, its applicability can be limited.

2.3 Origin of the speckle in the recovered image

To understand the origin of the speckle, we express the recovered object wave as follows

\hat{o} (x, y) = \sum_{k_{x}} \sum_{k_{y}} \hat{a} (k_{x}, k_{y}) p (x, y; k_{x}, k_{y})

= o (x, y) + \sum_{k_{x}} \sum_{k_{y}} e (k_{x}, k_{y}) p (x, y; k_{x}, k_{y})

where

\hat{o} (x, y)

is the recovered object wave and

e (k_{x}, k_{y}) : =

\hat{a} (k_{x}, k_{y}) - a (k_{x}, k_{y})

is the angular spectrum estimation (ASE) error for a given angle

(k_{x}, k_{y})

. Here, it is seen that the elements of the ASE error combine plane waves, and the combined waves are added in the recovered object wave in addition to the original object wave. Note that the combination of the waves generates a SF. That is, the speckle in the recovered object image, which is the intensity of the recovered wave, is made directly from the ASE error. The energy of the speckle is expected to become smaller when the ASE error decreases. With this finding, we come up with a new speckle suppression approach via reducing the ASE error rather than the ensemble averaging approach.

3. Phase conjugation method in TM-based image recovery

By using matrix and vector notations, [Eq. (2)] for the measurement model can be expressed as follows

y = T a + n

where

y \in ℂ^{M}

,

a \in ℂ^{N}

,

n \in ℂ^{M}

are the vector representations of

y (x, y)

,

a (k_{x}, k_{y})

,

n (x, y)

, and each column of

T \in ℂ^{M \times N}

is the vector representation of the

t (x, y; k_{x}, k_{y})

for a given

(k_{x}, k_{y})

, thus T is the concatenation of the vector representations of

t (x, y; k_{x}, k_{y})

for N number of

(k_{x}, k_{y})

angle pairs. Here, M is the number of measurements and N is the number of angular spectrum elements of an object.

Now, we review the estimation by PC in the previous TM-based recovery process [3–7]. PC uses solely the correlations between the columns of T and y for the estimation. The estimate by PC in [5,9] is written as

{\hat{a}}_{PC} = T^{*} y

= T^{*} T a + T^{*} n

where ()* denotes the conjugate transpose of a matrix. PC may work well when the columns of T are mutually orthogonal. In this case, the estimated angular spectrum becomes exactly the original angular spectrum,

{\hat{a}}_{i} = a_{i}, 1 \leq i \leq N

for

T^{*} T = I

; assuming noiseless situations. For correlated cases though, each element of the estimated angular spectrum is contributed not only from the angular spectrum element with the considered angle but also from those with the other angles whose SFs are correlated to that with the considered angle. Thus, erroneous estimation is made even in noiseless cases. For noisy cases, the estimated angular spectrum is corrupted further by the noise. Thus, the estimate by PC is expected to be erroneous in practice.

Note that turbid media do not provide orthogonal TMs for they have memory effects among the SFs of the input waves whose incident angles are not separated enough [17]. TMs may not have highly correlated columns for they have randomized internal structures; the degree of column correlations can still be significant. In our experiment, many columns of TM have correlations between adjacent columns around 0.5, for the normalized columns with M = 17442 and N = 20000. Considering a more underdetermined system, the correlations become severer. Thus for correct recovery, the correlations between the columns of TM cannot be ignored. This is why the process of PC recovery relying only on the correlations between columns of T and y is not desirable. Classical pseudo-inversion (PINV) approach can be considered for estimation. However, the inversion of the large TM requires relatively long computation time and the effect of measurement noise cannot be handled.

4. Sparse representation for speckle suppression in TM-based image recovery

SR has received great deal of interests recently for it enables correct estimation of a signal even for underdetermined measurement systems (M ≤N). Here, the measurement matrix surely has correlations which are not negligible among the columns of the measurement matrix. For successful applications of SR, there are two key conditions to be met: i) the compressibility of the signal and ii) the isometry of the measurement matrix. Here, compressibility means that the signal is well approximated with a small number of nonzero elements, say K where K = N. The isometry is the characteristic of the measurement matrix which preserves the energy of a signal well after the signal is transformed by the matrix.

Let us discuss the compressibility of signals and the isometry of TMs in TLI. i) Most natural object images are well approximated by only several terms in the Fourier domain [10]. We see that the basis signals in TLI are plane waves with different angles and the image is an angular spectrum in the Fourier domain. Thus, we expect most of considered images in TLI to be compressible. ii) Checking the isometry of a matrix is a non-deterministic polynomial-time (NP) hard problem. But, the Gaussian distributed matrices are proven to have an optimal isometry [10,18]. Through the random walk analysis, it was found that the SF in the transmission geometry is complex-valued Gaussian distributed provided that the number of elementary contributions is large [16]. This result is applied for the wave propagation through turbid media. Therefore, the TM is a complex-valued Gaussian and thus an effective measurement matrix. Now with the two conditions ‘i)’ and ‘ii)’ satisfied, we expect the benefit, a considerable reduction of the ASE error and hence the suppression of the speckles, in the TM-based recovery.

The K largest element approximation (K -LEA) of a signal is known to be optimal among all K element approximations in terms of the L1 norm and the L2 norm error senses [18], where the L1 norm and the L2 norm of an error e are defined as ${‖ e ‖}_{1} : = \sum_{i} | e_{i} |$ and ${‖ e ‖}_{2}^{} : = \sqrt{{\sum_{i} | e_{i} |}^{2}}$ respectively. The K -LEA can be made by an oracle estimator which is assumed to know all the coefficients of elements of the original signal exactly. It can be done by setting all elements to zero while keeping the K largest elements intact. The so-called sparse signal estimation (SSE), used in the SR framework, comes with a surprising result that it behaves like the oracle estimator; see [18] and references therein. In particular, the estimation error of SSE is bounded by that of the oracle estimator, as long as there is enough number of measurements. For a Gaussian TM, about $M = O (K \log (N / K))$ number of measurements shall be good enough [18]. This means that the SSE offers the perfect ASE in the TM-based recovery when the angular spectrum of object wave has K or less nonzero elements. For object waves with more than K nonzero angular spectrum elements, if the angular spectrum elements other than the K largest elements are insignificant, the ASE error of SSE is negligible. This is true for most natural object waves which are compressible. Besides, the effect of the noise on the error of the SSE is well bounded [18]; not amplified as it is in the case of PINV. Therefore, the SSE in SR frameworks makes it possible to lessen the ASE error considerably, hence to reduce the speckle in the recovered image significantly.

This SSE, an oracle-like estimation, can be done by solving the following L1 norm minimization problem [10,18]

{\hat{a}}_{SSE} = \underset{a}{\arg \min} {‖ a ‖}_{1} subject to y = T a .

The L1 norm minimization

{‖ a ‖}_{1}

promotes the estimate to be close to a compressible signal. There are many algorithmic approaches to solve the problem in [Eq. (8)]. We use an alternating direction method (ADM) [19] which is known to be efficient for the type of the problems represented by [Eq. (8)].

5. Results

We compare the angular spectrum estimation capabilities and the speckle suppression capabilities of three methods: two conventional methods (PC and PINV) and a SR based method (SSE). For fair comparison, additional information other than T and y is assumed not available.

Consider Fig. 2 where the estimated angular spectra are given for an object sample. It is found that the ASEs by the conventional methods, PC and PINV, are significantly erroneous. Except the dc terms, the angular spectra are all severely perturbed. On the contrary, most error terms in the estimated angular spectrum by SSE are reduced considerably. This error reduction is done while the features of the object image, which are the horizontal and the vertical patterns in our experiments (refer to Fig. 3), are all contained in the estimated angular spectrum, which is around the axes in the Fourier domain. The ASE capability of SSE is found to be superior to those of the other two methods. With this improved ASE performance, we expect to achieve significant speckle suppression in the recovered image.

Fig. 2 Estimated angular spectra using (a) PC, (b) PINV, and (c) SSE, respectively. Here, M = 4389, N = 20000 and K is less than 147. All angular spectra are represented in log scale for better visibility.

Download Full Size | PDF

Fig. 3 Recovered amplitude images averaged over one, three, five, and seven samples and cross sections of them averaged over seven samples using (a) PC, (b) PINV, and (c) SSE, respectively. Here, M = 4389, N = 20000 and K is less than 147. Scale bar: 10 μm. Constrast-to-noise ratios (CNRs) are calculated in the subsets (red arrow lines) of the cross sections.

Download Full Size | PDF

The speckle suppression capability of SSE is clear. We show in Fig. 3 the recovered images of the three methods when the numbers of averaged samples L are one, three, five, and seven and the cross sections of them for L = 7. PC and PINV have severe speckles in the recovered images. For a single realization, i.e., L = 1, the images do not match with the original image at all. As L increases, they become better, but they are still perturbed a lot by speckles. Besides, the two have speckles with non-uniform background densities. If the high resolution details of the object are placed around the corners of the view where the speckle densities in the recovered images are higher, they may not be interpreted easily by observers. In contrast, SSE is shown to suppress speckle significantly well. Even with a single realization, i.e., L = 1, it allows a modest image reconstruction. For L = 7, the speckle is successfully removed in the recovered image by SSE. There is little speckle and the features of the object are almost retrieved.

Due to the speckle suppression capability of SSE, the peak signal-to-noise ratio (PSNR) performance of the TM-based recovery is improved considerably. For example at L = 20, we see that the recovered image by SSE has PSNR around 60 dB, while those by the other two are around 50 dB (Fig. 4).

Fig. 4 PSNR performances of PC, PINV, and SSE.

Download Full Size | PDF

It is known that the turbidity in biological tissues poses a strong limitation for optical bio-imaging. We tested the speckle suppression capability of SSE in cell imaging through biological tissue. We performed the TM-based recovery for imaging of a microglia cell extracted from the brain of a rat through a thick (450 μm) skin tissue slice from the belly of the rat. Figure 5 shows that the quantitative phase image of the microglia cell is recovered successfully by SSE with significantly reduced speckle. On the contrary, the recovered images of the other two have considerable speckle. Hence we showcased the SR based SSE performs very well even for imaging through a thick biological tissue.

Fig. 5 Recovered quantitative phase images of a live cell averaged over one, three, five, and seven samples using (a) PC, (b) PINV, and (c) SSE, respectively. Here, M = 696, N = 5000. Color bar: phase in radians. Scale bar: 10 μm.

Download Full Size | PDF

6. Conclusion

In conclusion, we demonstrated a successful combination of the TM-based recovery with the sparse representation method. This approach enhanced the accuracy of image reconstruction and attenuated the effect of measurement noise. As a consequence, we could significantly reduce speckles in the recovered images and reconstruct high-contrast object images with much smaller number of raw images than that required for the conventional approaches. In this respect, our study will be beneficial to many of the recent interesting studies based on the transmission matrix. The reduced number of raw images will lead to speeding up the image acquisition and facilitate the observation of dynamics of the sample.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (Do-Yak Research Program, No. 2013-035295).

References and links

1. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics 6(5), 283–292 (2012). [CrossRef]

2. Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics 2(2), 110–115 (2008). [CrossRef] [PubMed]

3. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. 104(10), 100601 (2010). [CrossRef] [PubMed]

4. M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express 18(4), 3444–3455 (2010). [CrossRef] [PubMed]

5. Y. Choi, T. D. Yang, C. Fang-Yen, P. Kang, K. J. Lee, R. R. Dasari, M. S. Feld, and W. Choi, “Overcoming the diffraction limit using multiple light scattering in a highly disordered medium,” Phys. Rev. Lett. 107(2), 023902 (2011). [CrossRef] [PubMed]

6. S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nat Commun 1(6), 81 (2010). [CrossRef] [PubMed]

7. T. R. Hillman, T. Yamauchi, W. Choi, R. R. Dasari, M. S. Feld, Y. Park, and Z. Yaqoob, “Digital optical phase conjugation for delivering two-dimensional images through turbid media,” Sci Rep 3, 1909 (2013). [CrossRef] [PubMed]

8. I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, “Exploiting disorder for perfect focusing,” Nat. Photonics 4(5), 320–322 (2010). [CrossRef]

9. Y. Choi, M. Kim, C. Yoon, T. D. Yang, K. J. Lee, and W. Choi, “Synthetic aperture microscopy for high resolution imaging through a turbid medium,” Opt. Lett. 36(21), 4263–4265 (2011). [CrossRef] [PubMed]

10. A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” SIAM Rev. 51(1), 34–81 (2009). [CrossRef]

11. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. U.S.A. 109(26), E1679–E1687 (2012). [CrossRef] [PubMed]

12. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods 9(7), 721–723 (2012). [CrossRef] [PubMed]

13. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17(15), 13040–13049 (2009). [CrossRef] [PubMed]

14. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58(6), 1182–1195 (2007). [CrossRef] [PubMed]

15. T. Ikeda, G. Popescu, R. R. Dasari, and M. S. Feld, “Hilbert phase microscopy for investigating fast dynamics in transparent systems,” Opt. Lett. 30(10), 1165–1167 (2005). [CrossRef] [PubMed]

16. J. W. Goodman, “Some fundamental properties of speckle,” J. Opt. Soc. Am. 66(11), 1145–1150 (1976). [CrossRef]

17. I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett. 61(20), 2328–2331 (1988). [CrossRef] [PubMed]

18. E. J. Candès, Y. C. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant dictionaries,” Appl. Comput. Harmon. Anal. 31(1), 59–73 (2011). [CrossRef]

19. J. Yang and Y. Zhang, “Alternating direction algorithms for L1-problems in compressive sensing,” SIAM J. Sci. Comput. 33(1), 250–278 (2011). [CrossRef]

Speckle suppression via sparse representation for wide-field imaging through turbid media

Abstract

1. Introduction

2. System description

2.1 Experimental set-up

2.2 Image recovery using transmission matrices

2.3 Origin of the speckle in the recovered image

3. Phase conjugation method in TM-based image recovery

4. Sparse representation for speckle suppression in TM-based image recovery

5. Results

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Equations (8)

Optics Express